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ABSTRACT 

A mathematical model has been developed for the free 
radical bulk polymerization of methyl methacrylate initiated 
by 2-2' azobisisobutyronitrile. The dynamics of the system 
are described in terms of monomer conversion, zeroth and second 
moments of the molecular weight distribution. The kinetic 
parameters of the model are estimated using the experimental 
data of Balke (1972). The sensitivity of the output responses 
of the polymerization system relative to the kinetic parameters 
is investigated. 

The terminal control problem of obtaining a polymer product 
with specified conversion, number average molecular weight 
M and weight average molecular weight Me is solved by con- 
Sidering the polymerization temperature, and/or the initiator 
feed rate as control variables. The optimal profiles are 
calculated via the application of the Pontryagin Maximum Principle 
to the developed mathematical model. Two algorithms (discrete 
control method and shooting method) are implemented to solve 
the resulting two-point boundary value problem. A variety of 
products with desired molecular weight distributions are 
produced by applying the optimal control policies to the 
polymerization reactor. Furthermore, it is shown that it is 
possible to reduce the polymerization batch time by successive 
minimization of the objective function while keeping the 
same polymer quality. The extensive simulation results presented 
in this thesis clearly demonstrate the benefits that can be 
gained from the application of the optimal control theory to 


the industrial polymerization reactors. 
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CUAPICE Reel 


INTRODUCTION 


TerobyecCei ves 

The manufacturing of polymers occupies a significant 
portion of the chemical industry today. The polymer indus- 
try has grown tremendously in the past four decades as the 
uses of polymers have multiplied. Nevertheless, the pro- 
duction of a great number of polymeric products is even 
LOdaveanvart.. Resbasic reasons for athis are: 

CI)Minewareticuliveconcharacrerize and 
especially measure on-line the physical 
properties (e.g. molecular weight dis- 
tribution) of the polymer product 

(ii) The lack of good mathematical models 
for andistrial polymerization teactors, 
whiten makes ther application of /optimi— 
zation techniques to these systems ques- 
tionable. 

The present work was undertaken and motivated by a 
need to (a) develop a suitable mathematical model for the 
bulk polymerization of methylmethacrylate (MMA) and (b) 
Showstne feasible applicaci ons Ommoptimalacontro lato poly - 
merization reactors. Based on the above needs the goals 
of this work were defined as follows: 

(1) Obtain, a mathematical model for the free 


radical batch polymerization of methylmethacrylate. 
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(ii) Estimate the unknown kinetic parameters 
of the model from available experimental data. 
(iii) Apply the Pontryagin Minimum Principle 
to the reactor model and establish a systematic procedure 
LOrethestime optima licontrol of thesbateh process. 
(iv) Study the sensitivity of the optimal poli- 
cies with respect to the model parameters. 
The system chosen for this study was the free radical 
Do Weapelyuertzactongou MMA®initiated=pyea 2 9. 25) AzOD1Si-— 
SOputyLoNcuri le @CAIGEN) Catalyst,  *this polymerization 
was selected because of its commercial importance and the 
extensive experimental data available in the literature 
CIMEMeEDIyCiCa le eDLOper les andes Kinetic characteristics Of 
the system. 
This Enesice 1s Ofcanized Into: seven chapters... Ihe 
fi Gseycnapter states the objectives of this work and pre- 
sents a literature survey on the kinetics of the polymeri- 
Zoriouson MMAS and on the application: of sopeimal control 
fNeCOutes CO SChe ODLIMmMzZallon OL POlLymerI zation Systems. 
The second chapter deals with the modelling of the 
bulkepolymerization of MMA Tin a batch reactor. Using 
the general population balance equations, a mathematical 
model is developed that describes the dynamic behaviour 
of the system. The model consists of four non-linear 
ditterential equations and can predictymonomer, conversion 
and the first three moments of the molecular weight dis- 


tribution as a function of the polymerization temperature 
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and initiator concentration. The kinetic parameters of 
the model are estimated from experimental data. The 
accuracy and sensitivity of the model responses relative 
to the numerical values of the kinetic parameters is also 
studied. 

The third chapter discusses the application of the 
Pontryagin Minimum Principle to polymerization reactors. 

In the fourth chapter, two algorithms are developed 
to solve the minimum time problem. It is shown that it 
is possible to produce a polymer product with prespecified 
desired properties (conversion, number average molecular 
weight, weight average molecular weight) by controlling 
the polymerization temperature in the reactor. 

The fifth chapter presents the optimal temperature 
profiles computed for different desired molecular weight 
distributions. The sensitivity of the optimal temperature 
policy with respect to kinetic parameters is also investi- 
gated. 

pvt nessixthschapeer, bouhetne carelyeterecd ‘rate 
into the reactor and the polymerization temperature are 
considered simultaneously as control variables. An algo- 
rithm is developed to solve the minimum time problem and 
get desired physical properties for the polymer. Optimal 
temperature and, catalyst concentration profiles are pre- 
sented. 

The seventh chapter discusses the significance of 


this work and its possible extensions. 
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WS2 S8Literatune Sunvey aiKinetics Gof ttherPolymerization 
of MMA 


The single most important property of a polymer that 
determines its quality and final use is its molecular 
weight distribution (MWD). MWD of a polymer and monomer 
conversion ®int@atreactor ‘are tthe ‘results sofhtammumber lof 
elementary reactions occurring in polymerization and of 
their respective rates. In bulk free radical vinyl poly- 
merization, much work has been done to establish the 
Geact Hon mechanism vat = low .conversionsr(lesssthaatl0;con- 
version). Beyond this range, diffusion control of some 
reactions (e.g. termination reaction) coupled with experi- 
mental difficulties has severely hindered progress in 
kinetic modelling. 

Matheson and co-workers (1949) measured as a function 
of temperature the average lifetime of polymethylmethacry- 
late radicals in the photosensitized polymerization of 
the liquid momoner. By combining their experimental re- 
Sultsmwithne chose Of Schull zeandeb laschkemtt 942) ee@and Schulz 
and Harborth (1947), they obtained the rate constants for 
propagation, termination, and transfer to monomer for con- 
version up’ to 10%) "They asserted that) the accelerated 
polymerization rate) occurring in the later stages was due 
to a decrease in the termination rate constant. 

Nandi (1957), Ferington and Tobolsky (1957) pre- 
Sencedhrescul tes one thes bulk polymerizatrionvomMins AvThe re- 


ported experimental results were obtained with a wide 
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variety of initiators and under different isothermal con- 
dition seyeirommthase datal/ithey. calculated the rate constant 
for chainsitrans fersto: monomer! and ther ratio (ke 2/k,) 

where Ko and k, are the rate constants for propagation and 
termination reactions respectively. Their results were 

in good agreement with the results of Matheson et al. (1949). 

Hayden and Melville (1960), Tonoyan et al: (1973) 
measured the increase in temperature in an adiabatic poly- 
méerazationrreactor asi a;functionmofstime:y Theye reported 
that after 10% conversion, the polymerization rate and 
lifetime of the polymer chains increased due to the in- 
crease of the viscosity of the reaction mixture resulting, 
thus, in adecrease of the collision rate of the growing 
Padicalseetiney(foundioutithat,@beyondr407,econversion, 
the activation energy for the propagation reaction in- 
creased while the corresponding velocity coefficient de- 
creased, due to the monomer addition step becoming diffu- 
s#oneconérolled: 

Paul et al. (1973) observed a similar decrease in the 
propaea Gaoneratenconstant mand dstudiedcalso tthecing luence 
of thesstoichiometry of the initiator system (benzoyl per- 
oxyd@. Lauroyleperoxyde peN=Nedimethyl-p-toluidine) on the 
rate of polymerization. 

Pavlinets ‘and Lazar ‘(1973) investigated the polymeri- 
zation of MMA initiated by hydroperoxydes-S0, organo- 
complex isystems to high:degrees of conversion tat 10-45°C. 


They showed that the initiating systems used ensured high 
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polymerization rates even at room temperature with an acti- 
vation energy of 6.8 kcal/mole. 

Yokotavandyco-workers (1968) carried out rotating 
sector determinations of radicals lifetime for ten metha- 
erylates undergoing radical polymerization at 30°C. The 
rate constant of propagation k., and the rate constant of 
termination kL have been computed from these experimental 
results. However, there appears to exist a discrepancy be- 
tween their values and those calculated by other investi- 
SaAtors:, 

Balke (1972) investigated the bulk free radical poly- 
merization of MMA to high conversions. The experimental 
conditions of Balke were similar to those used in indus- 
bolalmreactOrc(nighvinitiator concentration and hich 
temperature). Using gel permeation chromotography (GPC) 
and a novel technique for the interpretation of the chro- 
matograms, he was able to follow the changes in the MWD 
as polymerization proceeded. He found that the rate of 
polymerization was proportional to the monomer concentration 
au bow conversion (less: than 20,)),.anGd ..oO.lowed a second 
order expression in monomer concentration after the onset 
Ofstnie zelrvetrect, Using Sawadase(1963)) equation co 
account for the gel erfect, Balke predicted the change of 
conversion with time. It should be noted that his model 
included two different expressions for the rate of poly- 
Merizeation, a LitSteorder Expression with respecc to mono- 


mer concentration before the onset of the gel effect, and 
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a second order expression after the onset of the gel effect. 
This introduced a discontinuity in the conversion profile 
plotted as a function of time. Apparently this model cannot 
be used for optimization studies. 

Mahabadi and Meyerhoff (1979) proposed a new model for 
estimation of the characteristic rate constant for primary 
radical termination, using the radical lifetime rate of poly- 
Merization and rate Or initiation. By applying this model 
to high conversion polymerization experimental data, they 
Showed the conversion dependence of the termination rate 
constant and the conversion dependence of the initiation rate. 

In the present work, a mathematical model for the free 
radical bulk polymerization of methylmethacrylate has been 
developed that accounts for the diffusion controlled termi- 
nation and propagation reactions. The kinetic parameters of 
the proposed model have been determined by fitting the experi- 
mental data of Balke (1972). Continuous differential equations 
which describe conversion, number average molecular weight, 
and weight average molecular weight are developed and solved 
numerically. 

Ivo leceratine scurvey: Optimal Control 

The Principle was first hypothesized by the Russian 
MaLtiemaercLan PONntuyagine ii) lo>O. ew inelooo, lt Wwase Lully 
proved by Boltyanskii and co-workers that the maximum 
principle was a necessary condition for optimality. 
Relatively little work has been done with respect to the 
optimal control of polymerization systems. 

Ray (1967) drew conclusions on optimum settings of 


temperature and catalyst concentration which minimize a 
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performance index described at final time only by conversion 
and polydispersity, in a series of tank reactors in which 
polymerization of styrene was taking place. He applied 

a peak-seeking method to compute the best temperature and 
the best catalyst concentration in each tank. Appropriate 
weighting factors were applied to each term of the objec- 
tive function in order to emphasize the importance of one 
variable over another. Hicks and co-workers (1969) formu- 
Patedsancimilarnsopvective function, and tried to optimize 
Loe sperrormance OL a polymerization reactor: (a) by con- 
trolling the polymerization temperature while keeping a 
Constant sinitietor concentration, (b) by using both tem- 
Deuatice anol tlator ‘CONnCenLrabions as control variables, 
The latter attempt proved unsuccessful due to computational 
problems. In both situations, the minimum time problem 
was not addressed. 

Osakada and Fan (1970) calculated sub-optimal temper- 
giucesand Catalyst feed rate policies an an attempt to 
obtain a desired molecular weight distribution. These 
near-optimal policies were represented by two time vary- 
Tee pOLwiomials the coetlicients of which were estimated 
Dyeeupactern search technique, coupled with a non-linear 
Saarcha technique. They found that the Suboptimal tempera- 
tiremand Catal vat tecdsrate policies were oscillatory. 

Greccitelli and Nicoletti (1973) cerived 4 modified 
form of the discrete maximum principle for systems with 


finite memory and presented a computational procedure to 
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find near-optimal policies in practical applications. They 
gave a numerical example relative to the computation of 

the optimal temperature profile in a batch reactor where 
two consecutive irreversible reactions were taking place. 

Sacks and co-workers (1972) applied the continuous 
maximum principle to minimize the reaction time for chain 
addition polymerization in batch reactors, while the nun- 
ber average molecular weight at the end of the reaction 
Was#hixedseulhesoptimal policyswas fournd®*to beva rising 
benlperature Sprofilemfor a tgiventeconstanteinitiator écon- 
centration: 

Shatkan and Gilman (1966) applied the method of 
variational calculus to search for an optimal temperature 
policy for batchwise thermally initiated bulk polymerization 
Slpseycene Welnathe searching *process® they *lirst .sétsthe 
final conversion at 95% and predetermined the number average 
molecular weight at some specific value. The optimal tem- 
perature so obtained increased with monomer conversion 
and reached an extremely high level (far above 200°C). 

Then they applied this temperature policy to the case where 
monomer conversion was nearly complete (99%) by following 
the optimal temperature trajectory until a predetermined 
upper limit of temperature was reached, then held the 
pemperaécuré@avethistlimit Seni ethesecndeoiethesreaction: 
Obviously the temperature policy for the entire reaction 
was not optimal. 


Show-An Chen and W.F. Jeng (1978) studied the minimum- 
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Gimerproblemeror the polymerization of styrene in a batch 
reactor by considering the polymerization temperature T(t) 
ands thetainitraleminttiacor concentration [I,] cyst GEG(ou Velehahwraeyih 
VaEvaoweses =  nelmekineric mode! ineltided the gel ertect. 
Calculations showed that the optimal temperature policy 
weacmoleNiticalm yepecter than the best Isothermal policy. 
Experimental verification of their theoretical findings 
revealed that there was a good agreement between experimen- 
tal@and@caleulaved=tinal conversion values. However, an 
appreciable deviation existed between experimental and 
calculated number average chain length. In 1980, Chen 
and Lin considered a two-stage polymerization process. 
PTResL vest ss ease or the=process, "initiated by a chemical 
initiator, was operated along the best isothermal policy 
computed to obtain predetermined number average molecular 
weight and monomer conversion at the end of the stage. 
The second stage was operated at a higher temperature at 
which thermal “initiation was important. They used the Mini- 
mum Principle to find the optimal temperature profile during 
the second stage, in order to minimize a performance index 
described by final time, conversion and number average 
molecular weleht “at the end of the™ second stage. Calculations 
showed that the optimal two-stage process was significantly 
béuser Gin’ terme or Lineal” time for a predetermined. quality 
Cfepolyeryrene)’ than’ che best isotliermal one-stage process. 
Masterson (Loy) mand Clough et al. (193) “applied the 


maximum principle to solve the minimum time problem for batch 
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polymerization of styrene. They obtained a product with 
desired conversion, desired number average molecular weight 
and desired weight average molecular weight by controlling 
temperature and initiator feed rate. The numerical solution 
has proven challenging. The search was successful only 
When Carried Out maErSe by a eradient approach, = then by a 
novel vertex technique. They found one temperature profile 
which conformed to industrial applied temperature profiles. 
iney neported ithan ay simian study would entail’ a cost of 

o/ O00 rine GPU “times 

Wu et al. (1980) gave a shrewd graphical solution to 
the Minimum time problem of styrene polymerization. Their 
performance index depended only on conversion and number 
average molecular weight at the end of the reaction. Thus, 
by applying the maximum principle, they were able to obtain 
analytical expressions for the costate variables and derived 
aQudiiserential equation descoribing: the optimal variation of 
temperature with time. Their theoretical predictions were 
in agreement with experimental measurements of conversion 
and number average molecular weight. 

Ine this’ work ia general alsorithm which Solves ‘the 
minimum time problem for obtaining a polymer with desired 
final conversion, final number average molecular weight 
and weight average molecular weight, is developed. Optimal 
Eenperatire polsciessand inleragor feed rate programs “are 
computed for obtaining a polymer with desired molecular 


weight distribution. A pilot scale batch polymerization 
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reactor has been built by the author to verify the optimal 
policies as part of another project. In this equipment, 

the temperature of the polymerization reactor is adjusted 

aS per predetermined optimal policies by supervisory compu- 

ter contuol MeL seunemtemperatire (loopyunder tanaloguase 1D icon—- 
trol, Conversion, number average molecular weight and 

weight average molecular weight can be measured on and off-line 
and checked against simulation results. A detailed description 


of the experimental apparatus is shown in Appendix E. 
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CHAPTER 21 


DEVELOPMENT OF POLYMERIZATION KINETICS FOR MMA 


II.1 Description of the Reaction Mechanism 
Polymethylmethacrylate (PMMA) can be obtained via 


the free radical polymerization of MMA monomer. 


CH CH 
3 3 
we | 
n CH = C - CH, — C-— 
mS ae 


(MMA) ae 


A peneraile description, ofthe reactionswhich are taking 
place durmnertne free radical polymerization of vinyl 


MOnOMerCS@imLliatede by ‘a free radicals catalyst iS"%as) follows: 


Catalyst decomposition 
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Termination by disproportionation 


tose eb 4h ’ (x,y) E[l,~) X[l,@) 


Termination by combination 


Re + RY Kte ee, to | ea Ste) Sanh) 
where Ifidenotes! the ainitwdator: molecules, Ry» Lhevya ni tiaitor 
Radwes sive nid R, , uhe acta viatvedi anivwelaitower ade calics lv 
denotes the monomer, Re represents the live radicals of 
Giaiieeenvthex eh is shel initiator eliaciency factor, Ka. 
k,, ky Ke, ka and ee Bree Gnenreacelconsteniss Don init iAcor 
Gissociation,;-initiation, propagation, transfer to monomer, 
termination by disproportionation and termination by com- 
bination reactions respectively. 

Dicainitvaror Wsed in this studyois 4 7 azopisiso— 
Diiewromliriie (ALBN) which is very common in, industrial 
Be-ClCrs pics decomposition rate, tollowesrirct-order kinetics: 

Using the above proposed kinetic scheme, equations 
GecCripineeLNemcOnversion Of monomer y Phne concentration 
Ommintiiaror.. and the molecular welent distribution, (twp) 
imap aLchyreactom ere derived, 10 Simplity the machemacical 
description of the system, the following assumptions are 
made: 

Gi) All reactions: are irreversible 


(ii) Reaction rate constants are independant of 


chain length. 
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(iii) Transfer to monomer is negligible. 
(iv) etermination reactionsie by dispropomtl onatuion: only. 
(v) Reactor contents are perfectly mixed and there 
ere eno stemperactire pradients in the reactor. 
These assumptions are common and well documented in the 


modelling offehe free radical polymerization of MMA. 


iiez, sDeriveatvonvor the Kinetic Equations 


The catalyst decomposition reaction follows first- 
ordersikineticsemd Ihe rate pofidisappearancer of catalyst is 


described by the following differential equation: 


(ee eee CEN) 


wIeweunl his the (nttTator concentration. (& 1s the time, 
Ka is the rate constant of catalyst decomposition, ne the 
Catealyss feedtrate in” the reactor during the polymerization 
and V is the volume of the reactor. Any variation of the 
volume of the reacting mixture due to density changes and/or 
Hee tasCr Pauaition is assumed to be negilivable. 

Pertinent Steps in the derivation Of. the equation 


describing the rate of disappearance Of monomer are outlined 


in Appendix A. The monomer conversion X is defined as: 
[M] - [M)] rAbys 
X a pe aaa = 


where {[M] is the monomer concentration in the reactor during 
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the polymerization, and [Mo] is the monomer concentration 
in the reactor at time t = 0. The change of the monomer 
conversion is described by the following differential 


equation: 


aaa Ky (1-X) CLL) 


where £ is the initiator efficiency, k_. is the propagation 


P 
Faces constant and Kk, US =the Lermination 2+ate constant. 

The MWD of polymers can be expressed in terms of its 
Momenes when the MWD FS unimodal, as it 2S the case for the 
DOlymerization OG MMA in’ the=renee Of interes: (Balke, "L972 je 
The moments of the live and dead polymer distributions are 
defined as follows: 
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where M1 ies: Glew moment of the Live polymer distribution 
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and [Re] fie the concentration Oo. Vive radicals er chein 


ben mies edil tle GeaccoL. 
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where uy is the fae moment of the dead polymer distribution, 


and Peg ee the concentration of dead polymer of chain 


ch 
lengthex in the reactor... The zero” moment of the live 
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polymer distribution, Ag» represents the total concentration 
of live polymer in the reactor. The zeroth moment of the 
dead polymer distribution, Ug» represents the total con- 
Gentration of dead@polymer in the “veactor. Theksum 

(uy taq) represents the portion of monomer which has reacted 


and=sis equal to: 


X (Fs A0)) 


The differential equations describing the variation of the 
moments of the dead polymer distribution are derived in 
Appendix A. The rate of change of the zeroh and second 


moment of the dead polymer distribution can be expressed as: 


a= 2 Ekg [T) (ity) 
2 
du k 
ee (II.8) 
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Diysical properties Ont epoiymer PGsuch js melingsas- 
COstty sechenieal résistance Silexip1lt ty eis tr engi hives) 
do depend upon the number average molecular weight, M_, and 
the weight average molecular weight ore of a polymer. ML and 
M are directly related to the moments of the MWD. 
ie al 
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Uo =f No 


where MW is the molecular weight of methylmethacrylate 
monomer. Since hie is negligible compared to UR Oy << Wy); 
equations CLE. 3)eand s(11-10)scan) be simplitied by neelecting 
Ehescontri bution of ho» hyo and ho moments to the calculation 


of Mo and Hc Therefore, 


{ = We — = MW Citta”) 
* 20 0 
uv 
PES Spall 9 
Mo = MW = MW x [Np] COE) 


EbemPpOlvdaispersity of a polymer CPD)a6 defined as the ratio 


GMa), ang as ene indicartion of the breadth “of the MWD 


PD: = = ———— (eas) 


In bulk polymerization of methylmethacrylate, dramatic 
physical changes take place during the course of the reaction. 
As polymer concentration increases, a point is reached 
where appreciable chain entanglements occur, and eventually 
a glassy state may result. These physical changes often 
Mave a signiticant €rrecton Dou rate OL polymerization 
and molecular weight development, and any attempt of model- 
ling such reactions must properly account for these phenomena. 


Experimental data show evidence of an autoacceleration 
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of the rate of polymerization at conversions above 10 - 20%. 
pecause OL thes increasing viscosity of the mixture, the rate 
of termination becomes diffusion controlled. Polymer 
chains begin to entangle causing a dramatic reduction in 
Ladreca | eChaimemobiiaty, piving. ay drop in kL. Because of 
ENeuNUCe her RCOnGceneLacion.OL free) radicals anethe Geacton, 
Piewravce, OL. polymerization increases. lt  ispinteresting 
LOunOLe sLiareevenerelatively slow propagationsreaecriens 
involving the small monomer molecules become diffusion 
Concrolled: | nev propagation rate constane KD drops 
dramatically. That means that no more reactions are taking 
piace in the reactor, therefore a limiting conversion below 


100% may be observed. 


Pit mee CoMme ler bel imaticonm andbransnereneocy sibivity 

The parameters which appear in the kinetic equations 
(ese wel lesen lao ear Ge the gin ite COnmeriaClency, she 
Biewrace, Constant, Of catalyst dissociation: Kae and the 
Vaetouct LoCeScuase Of Uhempropagation race cons lane ever 
LNewLermineclOon rete Constant Ky Ca) kD. 

THheglMeLlaconvetricuency, Bias Lenhescemis the Lracmi0n 
CMe nececOta le numbDermOL mini tlotOmeuad 1 Cals atinadt gat Gxaclua ily 
“ised .in the tormation of polymergpchains, fis ,an,empinical 
GAaCtOr with vVelués in the ranges 5 -sl0s lhe selriciency 
factom usually talls withbeincressinegsreacti on wtemperactunes 
and with conversions of monomer to polymer. The changes 


however are relatively .small and it is reasonable toassume 
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that £ is independent of temperature and monomer conversion 
PorspracLi calecarcilations, Iinethistwork,. thesvalne of. f 
has been chosen as 0.6, figure commonly used in literature. 
The decomposition rate constant of AIBN, ka has been 
Gdecermaned byamany= worker's .~ Da tava htersiba lWel (107205 
Or Dratscoll, and’ Dickson: (1968 a0 and@iitor {196 9998 were’ cornel— 
ated with an Arrhenius type plot giving the following equa- 
tion which was used in this work: 
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d 
where T is the absolute temperature in the reactor (°K). 
As explained above, the termination rate constant, 

kL, and the propagation rate constant, og change with 
Cconverorvon.§ / Asse fresu eof this-theltatdo (Ky = ke/k.) 


Wiel eso Wary with Conversion. ©'To? express. thst variation 


qiantirative ly the following tequationis wsed 


Ky = Kog CLS Me eee ay) GRE ty) 
where Koo is the value of Ky at conversion X = 0. g(X,T) 
tray be called» the "els effect function, because it accounts 
for the conversion dependence of k, and Kk. due to the diffu- 
Sion control, of termination and propagacion zeactions. 

Koo can be determined by measuring the rate of poly- 
merization at low conversions. Data were correlated in an 


Arrhenius type plot by Balke (19/2).> which) gave tae follow— 
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ing expression for Ko: 


Koy = 100.4 x EXP (-2960/T), (1.mole’'s™*) (11.16) 
Because of the gel effect the termination rate constant 
decreases silhisszresulits in “an increase of the rate of poly- 
meri 2a blo miyethen conversion mance «oh 20 — 6077 » Beyond 
this conversion range, even the propagation reaction could 
become diftustons controlled, which inwturn could cause a 
Cuamatd cuGecrease onerhe, polymerization rate... To,describe 
this unusual behavior of the polymerization rate over the 
whole monomer conversion range, the gel effect function, 
e(X,T), must first monotonically increase up to some monomer 
conversion value (60 - 80%, depending on the polymerization 
temperature), then decrease up to the final limiting con- 
VelLelOnevalucape 1 OveCCount fOr sthisetypesor behavior, gan 
exponential relationship of g on X and T was assumed. Friis 


ana Hamielec (19/5) have used a similar expression. 
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Pte = EXP (AX POBKS ac) nia 
with A = A,/T + A» 
B = BL /T + Bo Giga SC») 
C = C,/T + Cy 
Ay> Bi» C1, Ay) Bo, Co, are unknown parameters estimated by 
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fitting our model to the experimental data (conversion, 
Mo Mi of Balke obtained at different temperatures. 
To obtain estimates of these parameters (Table II.1), a 


Finite difference Levenberg Marquardt routine from the 


IMS Library was used: 


LRapvera tert 


=96900 = 50. 


Maci Gare Nac £O De Laken in tne. estimation, of A, &, 


C parameters and especially in the choice of the initial 
guesses of Ay» Bi» Ci. Ay, Bo, Co, because the system res- 
ponses are sensitive to the numerical values of these para- 
Mmevers, we Lioetetsls -dtscissed Later in his) section. 

The results of the parameter estimation are shown in 
Prices ive, CO LE 6. sin these fapurecee Calcuimeec results 
for monomer conversion, MA and ag are compared with the 
experimental values of Balke obtained at three ditfterent 
POLyMNer uzation temperatures.” “Ihe solad Lines represent the 
calculated model response and the small squares represent 
Chee ceperimehtal pOlntc, Om bathe Gly 2)re Bl urcaumDcescen 
Bhat the model ts describing accurauemy tne Vvaldatlons or 
conversion with time. The slight discrepancy observed 
Iieeteeine the molecular welgic wdasiE EpUi lor (MO oe could 


be due fo possible errors in e€xperimencal measurements. 
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Bicune sl seuspOwomcrhe.variation) of athe funcéeion je OT) with 
conversion at two different temperatures (50°C and 90°C). 
Although the mathematical model chosen in this work 
can describe satisfactorily the polymerization of methyl- 
methacrylate, yet our knowledge of the kinetic parameters 
of this model may be imperfect. This gives rise to the 
important problem of parameter sensitivity, which is de- 
Pipedmeasetieeecrrects Ofeuncertamties im theprabe coerticients 
on the calculated output responses, namely, the effect of 
pamaneater uncertainties ein Ka; wy Ky, on the calculated 
Paty Ree ay Ug» Ho: In. addiiwonto.taessensitivitysorethe sys- 
tem responses to the above parameters, we are also interes- 
ted) in knowing the effects of the initial catalyst. con- 
eentration [I] and polymerization temperature IT on ythe 
output variables. 
The sensitivity coefficient for the parameter, Pas 
and the output, Zs is cebineds as) tien 7 Sly At io ecer ay — 


ative of Z. heijela ey vel Noteyerenes (are Ds 


ie al 
eae Sap CLES) 


SenGitigity woekricients indicate the magnitude and the 
direction of change of the response Z due .to perturbations 
inegche values Of the parameters g-soince stbe model. sic.a 
Set, Of non Lineam Obdinany, dliler emia) Seca one sacice sen 
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cation of the well documented sensitivity analysis to the 
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model equations (Atherton et al., 1975; Beck and Arnold, 
PU) ery relvis. a system of *ditterential equations tor tne 
Sensrervyityecoeeurcients, Calledesensitivity -eduarion 
(see Appendix "8B “for “further details). ©°The variation of 
ier sens rorvyi ty =coctrrcrtente with respect co "time Ls ob- 
tained by integrating numerically the sensitivity equa- 
PPOs she resets mOlmLnie Intvesrationeot tiemsensitiyi Ly 
equations for ka PP Kyo» Ay, By, Ci; Ay, Bo, Cy are summ- 
Paes Ceti en te rcs mileso rand fermen order ro show erie 
relative influence of the parameters on the output varia- 
bles, normalized sensitivity coefficients, as defined by 


equation (11.20) have been plotted. 
ety =e So iuigae0) 


AepoOstt Ve Isensitive cOeriicient indicates that a posieive 
Vr teaenror OL Lie COLrrespondine PDaramecer Causes (an in- 
Creasewin thie output variable. A negative sensitivity 
ecgetricr1ent indicates: thae a positive variation of the 
parameter results in a decrease in the output variables. 
ite LS asec thatethaioutpue variables) iareimosh-cens peive 

to Ay» Bi» Ay, Bo parameters. 

An increase in Kg Orin ee dOeCeenOlcauscranyc12 = 
piercer perturpation On sLhesinitiatom concentration angson 
the seen! moment, however it causes a slight increase 
in monomer conversion, and a slight decrease in second 


moment at the end of the reaction (hicures 11-6 sand 11.9). 
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The effect of perturbations in the reaction temper- 
ature is shown in Figure II.10. An increase in temperature 
results in an increase in monomer conversion and zeron 
moment, a small decrease in initiator concentration, and 
a decrease in second moment at the end of reaction. 

The effect of perturbations in the initial concen- 
Eeetion Om Catalyst [Ig] is shown sin’ Figure (ali eee he 
output variables are much less sensitive to [Ip] thany.co 
T. However, the profiles of respective sensitivity coeffi- 
cients have the same shape. An increase in [Ip] causes 

2 
an increase in monomer conversion and woe moment, and 
a decrease in the second moment at the end of the reaction. 
Tmesquantitative interpretation ofwmigures! 11.8) te (11. ia 
Shows how accurate the estimation of the kinetic parameters 
should be. AS an example, consider the effect of a pertur- 
bation in B, on the conversion X&%. “Insiicguresl1.8, 1 1s seen 


that the dimensionless sensitivity coefficient Vag can be 


equaieup to 25, 
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Figure II.4. 
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Number Average Molecular Weight and Weight Average 
Molecular Weight versus Time._ Isothermal Run at 
50°C. (Io = 0.01699 gmole. 271) 

____ Mathematical Model 

© @ Experimental Points (Balke - 1972). 
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Molecular Weight versus Time, Isothermal Run at 
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____ Mathematical Model 

© cw Experimental Points (Balke - 1972). 
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Molecular Weight versus Time. Isothermal Run at 
90°c. (Io = 0.01699 mole. 27!) 
____ Mathematical Model 
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CA PLE hae 
APPLICATION OF THE OPTIMAL CONTROL THEORY 


TO THE POLYMERIZATION REACTOR 


Pie eeroml aclonsOrme le (Optimal Control Problem 


The theory of time optimal control has arisen from the 
need to improve the performance of dynamic systems. More 
specifically, given a process described by a mathematical 
Model lt is required to find the admissible inputs (Cor 
control variables) which generate the desired output and 
Which mineso doing.) Minimize a choSenwcocterunetional, 
Stated mathematically, problems such as these belong to 
the calculus of variations. However, classical variational 
theory can not readily handle the hard physical constraints 
usually imposed on state variables and control variables. 
Mise otericuley led Pontryagin (1956)" co rirse conjunceure 
his celebrated "Minimum Principle", and then together with 
Boltyanskii and Gamkrelidze (1962) Mtoe prOviaespLoon Ore ate 
The Minimum Principle gives a solution to optimal control 
problems which cannot be solved by the Class teal) Vartacionead 
LHeOLry. 

Consider any controlled process which can be described 


Dymeesy ocean. Ore ordinary differential equations: 
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where Y is the state vector of dimension n 
Y, is the initial value of the state vector at time to 
Ge elcetnes time 


u is the control vector of dimension m 
Given the initial value for Yor equacion. Gil, em trneisolieion 
for equation @It-1) is uniquely defined, once the control 
VeCHOr UW is Chosen. 


Consider the performance functional J: 
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where to and te denotes thevinitialvandstunal times G andl 
are scalar functionals. It is required to find an admissible 
control u*(t) which causes the general non-linear system 
CrEle i) with gnitial condition, (Lil. 2). te. follow an admissi- 
ble trajectory Y*(t) that minimizes the performance functional 
ee Heme SICH. 4a) Contr ou. u(t )e is: called the optimal control 
Busi) scoce Optimal trajectory, 
ties chassicalecalculus of varigtlonss | eadcato the tollow= 
Ties SOLU tion: 


Detine ihe Hamaltoniian H (scalar functional) as: 
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Satisfy the canonical equations, 
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The final boundary values for costate vector P depend upon 
the particular foun of G(Y(t-), te) and upon the final con- 
C@erTors TUT tie wotatrer vector Y(te). Titaiseiice 

(Py "1 the=costate vector Y(t) is free, then 


P(te) must sadlisty ene following ‘equation, 


3 G(Y¥(te), te) 
© = "OGD Tee Gisiis +6) 


P(t 
(ii) If some of the components of ¥(t¢) are specified 


aomene £t inal "times “ire 


a. 2, 


¥;(t¢) = Ye er eee?) sg 
Gis) 
ea ee is free EO) ee n 
then P(te) MOS taSabisny Che Tollowinesecondlelons.; 
P. (te) is free pee pet al ee te ae the 
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A fecessary condition ror optimality, derived from the 
Qlasstcabecarculus of Variations is that) the radient) of 
GhemHami | LOntaneH relative COs tice CONCLOLMVeCLOm Uusmus 
pe zero alone the optimal, trajectory: 
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40 
Asheuristicaproor Of conditions (211.6), (T11l-8) and 
(III.9) is given in Appendix C. Equation (III.9) is valid 
only if the control u is not constrained. However, in 
most control problems, restrictions have to be imposed on 
thei controlyuy (eng .4) u, Gt iia < M;, Joma \222.%% m)? 
because of economical and physical constraints. Then the 
above relationship (III.9) cannot be used because it breaks 
down on the boundary. This relationship is replaced by 
the following inequality: 


Ging ultd Pans HeseH (SF eu, Pigew) CEH 0) 


Poneal admissible: controls u(t), and all ate tel, 

Whemee Venice Prerepresents theeoptamal trajectories for 

leeUle re eeuetoniinimiZed by Sselectinesthe admissible controls 
which give the smallest value for the Hamiltonian Equation 
(III.10) and which is the basic statement of the Pontryagin 
Minami rincaples Solution to the optimal control problem 
can be obtained by solving a two-point boundary value problem 
Clon meaOme ie wditcerential equations. (hluel) sand iio). 
WLC ree iee lal values LOL Y(t 0% Cetinal values sor (te): 
and (n-r) final values for P(t,), while minimizing the 
Hamiltonian H. It must be noted that (a) it may not be 
ROW aavence, that, an Optimal” CONtCrOl exists and "(bD) 

evel 1 al OpLimal CONtrod® exists, #lt=may noc be ullidue: 
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reactor can be described by a knowledge of initiator con- 
centration fiat monomer conversion X, ae moment of the 
MWD Ho: and second moment of the MWD Uo: These state 
Variables, obey the adirlerential, equations (11.1), (T1-3); 
Chin /)- = (blo ymderived in Chapter ll.) Because the state 
VarLaprece| |i mx. Ug: Ho are of a different order of magni- 


tude ({L] = 1072 


se Olde he eke uy #10"), the numerical 
integration of the differential equations has been proved 
challenging and was eased by the introduction of the follow- 


ing state variables: 


Dimensionless catalyst concentration: Yy = Tne 
0 
: : . 38 
Dimensionless monomer conversion: Yj = ¥— 
d 
th a) 
Dimensionless zero moment : Y, = — 
Od 
me 
Dimensionless second moment: Yy = — 
aod 


where Xa, Vog? Hag are the desired tivalucs: com tx) HO and 
Uy a Gena itbime te. In terms of the dimensionless state 
Viamuake smequations:) (Bbal)ipaCi irs 2, Cie, DrlancirGit Seiten 
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dY, - 
qpeee char sate1at VITo] 
dY, (lene enka) 
Ace Ma CES ae lol Ko XK (Geemalay) 
dy 

3 1 
aa- = — 2£k, [I,] Y 
dt og ay. S02 a4 
dy 
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The initial values for the dimensionless state variables 


are given as 
(ayaa ce, to. Sto fhtos) (LT UAEZ) 


The choice of a mathematical performance functional is 
a highly subjective matter, as the choice of one engineer 
needs not be the choice of another. The objective functional 
is based on product quality and economical needs. In this 
particular problem of polymerization it is, desired, to 
produce the "best polymer" in the minimum batch time. 
The mathematical definition of the term "best polymer" 
depends on the final application Of fthe polymeric=producr, 
hence it can widely vary. The useful mechanical and phy- 
sical properties of a polymer are strongly related to its 
MWD. In most cases, there is a molecular weight range in 


which a polymeric product will have desired physical and 
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mechanical characteristics. Thus the term "best polymer" 
will mean a polymer with specified MA and Me, and conse- 
quently will possess certain desired physical and mechanical 
properties. Prespecified conversion, number average molecular 
weight and weight average molecular weight for the final 
product determine uniquely the zero, first and second morents 
of the MWD at the end of the reaction. Accordingly the objective 
function to be minimized may be expressed as: 

aes eR hs Wes (eetee) aed) s 


2 Gone 15 


oe (M7 Gor) a) 
This is a Mayer problem. 

The final time ty is fixed in the above objective function. 
By successively minimizing this objective function for smaller 
values for te, the minimal batch time problem can be solved. 
Comparing ecuations (Ilil.3) and (lii.13), we cbtaan othe 
following expressions tor L and G: 

Ore ce) = (iy (he) 1) ru oe 


pa 


<P (Y, (te) = 41) @E1T 14) 


Tia), Bee (Git) ee = tide 


Because the final conditions on Y have been included in 
equation (III.14), ¥(te) will be free, and therefore 
P(t,) will be given by equation (111.6). 


Tas to be minimized by finding the optimal policy for the 
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control u(t). In this study the two control variables are 
the reaction temperature T and the catalyst feed rate Ee 
It is assumed that Fo and T can be manipulated and fixed 
at our convenience, although they could be bound by physical 
Constraints. 

The expression for the Hamiltonian can be derived from 
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Heisenote an exDlL icles nCi Lone O canLumles 
THererore, as ils proven in Appendix C,H must remain icon= 
Stumeraloue the optimal «trajectory if sthevequation se (111.9) 
Lemcombes satist ied. 
The costate differential equations are derived from 


equations (III.5) and CU eae Ge 
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pier! #2 18, 

Po (te) = 2(Y5(t,) = 1) 

P, (te) = 2(Y, (te) =e) (GOL 18s))) 
P, (te) = 2(Y, (te) - 1) 


inemevccenswOLedirterentials equations (1iilal)) and) (lit a7) 
have to be integrated simultaneously with the initial and 
Pine eCOnC LCOS sClW le 2 )yeand (ill als )28) Theretone: ay two- 
point boundary value problem has to be solved, along with 
the minimization of the Hamiltonian H to obtain the optimal 
Contzol poligies. 
In this work three separate problems are considered: 

(GDR ODealneaepolymeLuwitnedecared 6X. Mek oy ae cBe inchs 
Minimum patch time by calculating the optimal temperature 
Demi Cyeticne Leactor. 

(it) Obtain a polymer with desired x; Mee in the 
minimum batch time by calculating the optimal policy for 
Cauauyst teed rate. 

(iii) Obtain a polymer with desired X, M., M, in the 
minimum batch time by calculating the optimal policies for 
polymerization temperature and catalyst feed rate. 

Solution for problem (i) is developed in Chapters IV 


and V. Solutions for problems (ii) and (iii) are developed 


ins Chapter V1. 
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CHAPTER IV 


SOLUTION OF “THE IPBVePROBLEM 


iV ee lneroductuon 

Te this chapter, we solve the TPBV optimal control 
problem by assuming only one control variable, the poly- 
merization temperature I in,the reactor: (‘Noi catalyst 
is added to the reactor during the polymerization, 
and the initial concentration of catalyst [Iy] 
toetaxede co a value or 0-5), in welghe.. ‘Two aleornithms that 
compute the optimal temperature profile in the batch reactor 
are developed and their results are compared. The first one 
Tse bases onethemaiscretizartion of the total’ polymerization 
time in N equal intervals. The control variable (temperature) 
remainseconstant in each time interval, and the Hamiltonian 
is’minimized by a first-order gradient method. For the sake 
OeEDUe vi mys icnis aleori cium wi lie besrerert cd lopacsenesiiccretre 
control method (DCM) algorithm. The second algorithn, 
which solves the two-point boundary value problem by a 
shooting method will be referred to as shooting method 


CoM)aleori chm, 


IV.2 Discrete Control Method Algorithm 


The polymerization of MMA in a batch reactor can be 


described by the following set of differential equations: 
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WLED thes initra lecondutyvonse: 


MAG) HS Ts (IV.2) 


The Objective function to belminimized is defined as: 


TNE GED) <1 AS ORGEAY = TY es ats) = TY 
ne 
+ We dt Give» 
=O 


ACCOLC ING StOBENemDEcCVIOuUS anaikyscis in Chapter Tlie the 
NeceSsaby, cOMdUL1Oons#Lor thes optimadaty of the costate 


trajectory can be expressed as: 


be je 

Le (IV.4) 
Ay LAOAGE A he ee) 

de) = a (IV.5) 


Suppose that a temperature history TS). t SUC tel is 
known and used to solve the differential equations (1V.1) and 
(IVe4) eeesoe that the Sstate-costate trajectories yi (t) and 
p(t) ea itier yethiemboundaryacondi cl Orca iVEZ mand (lVeo)m | Li 


this temperature history aie aa also satisfies equation ive o)G 


palaces pie 7 ie a) e= 0 (IV .6) 
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for alt (| ren tel, then Ta() satisfies all necessary 
conditions for being an optimal control. 
However, suppose that equation (IV.6) is not satisfied, 


then for a change 6 tT+(t) around the control T(t), 
eT (eo), 2a De ee PO (IV.7) 


the first variation of the functional J may be expressed as: 


(See equation C.12, Appendix C) 


Ue A 
r abst \ 4 i 
é J -| Ga chp ade (IV.8) 
tC 


O 
Hence, if the change in T+ is selected using a steepest 


gradient method 
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Fy] 


: ; i ; 
Suite? Sat (4) with a? > 0 (IV.9) 


Then 
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Because the integrand is positavyertor ali Peta tel, 

* j= 0. For a sufficiently small 6 J (ise., a Suri icrven tisy 
small oy this implies that a(t (t)) < I(T (t)), since 
the increment in J is governed esseptialive bye these nse 
variation 6 J. Hence, if the Contiow tht) is updated 


according to equation (IV.9), T(t) will eventually converge 
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touche optimal control T*(t), if the optimal control exists. 
We shall now outline the steps required to determine the 

optimal control T*(t) using the discrete control method 

algorithm (V. Gourishankar, 1980). 

Deepal: selectea discrete approximation to the temperature 

HESscory T(t), to oe ee te, by subdividing the interval 

[t., tel iNncOomNWesultably chosen subintervals of equal dura- 


tion. Consider the control T’ as being piecewise constant 


Guning seacn Of Enese subintervals, that 1s: 
i fi al = 
Ie) AN OS Sy eae a Bll kel S25 mca OYE TE) 


and store the value of Ge). 

Stcpes- 2 Us. e tht), integrate the state equation GLY» 

from to to te with the initial condition (1V.2Z) and store 

the resulting state trajectory ECS). 

Step 3: Obtain Ba (te) by substituting Miele) in equation 
Clive yuausingeathis value of (ee): the computed values 
Geethe state trajectory fey, and, thescontrolehistonry Tit), 


integrate backward equation (IV.4) from t, to t); compute 


; sh 
L a) ee (: | 
g (ty, 9 lon 


in each interval eases ty] and save this value. 
Step 4: If the stopping criterion is satisfied, terminate 


; vi : 
the iterative procedure; in thie casera GpaeGo eres rediLLod 
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optimal control. Otherwise generate a new piecewise constant 


eontrol history 
Ta ee ee) orate Mma, 2 Ne (TV TZ) 


and restart the procedure from Step 2. 
This algorithm is summarized in Figure IV.1. 

ite themabove aloor tim, Che sindl times faxed. 81 ne 
minimum time is found by solving the same problem for diff- 
erent final times and selecting the smallest final time 
for which convergence is possible. 

inthe. DCMealeond chm: TH(t) ce no tecOnscrained. 
theretore, the necessary condition for optimality relative 


POmLnesminimuzacion of H is 
o*(t) 3.) ian = 0 ‘gael kil Picaltes te] Gives.) 


Consequently, the above condition should be the stopping 
criterion for terminating the DCM algorithm. However, 
instead of using the above criterion, we have chosen the 


following stopping condition: 
eee yor (IV.14) 
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was found that when condition (IV.14) was satisfied, con- 
dition (1V.13) was also approximately satisfied. To verify 


that, the following norm for the gradient has been defined: 


Me th Lend 
Sy hh Gay) ORAS) 
k=1 


As the above norm decreases and approaches zero, the temper- 
ature profile computed by the DCM algorithm approaches the 
optimal temperature policy. 

The integration of the differential equatronss (ive) 
and (IV.4) has been carried out using a Runge-Kunta routine 
of Shampine and Allen (1973). The routine estimates the 
becCaiwerror and adjusts the integration step size depend- 
itebyson then localserror. The initial integration step 
size has been selected as 5 sec. Stiffness in the differ- 
ential equations was not encountered, because of the use 
of dimensionless variables. 

The gradient o (ty) Can be  calculaledeanalytically 
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The following numerical values were used for E, and E 


d 2 


(Balke, 1972): 


Eq = 32387. cal. mole™!; E, = 5882. cal. mole” 


ut 


The choice of the length of the time subinterval 
[t,_,> t,] is important. Figures IV.2 and IV.3 show the 
effect of the length of different subintervals on the 
optimal temperature profile. A longer time interval 
(25 min, 10 min) does not change the general shaper or the 
optimal profile, but produces an obvious dis Conta nussta et a 
PCs. shorter time intervals (5s, 1 Hin’ 2 Min. pman ) 
produce almost a continuous temperature ProOrillen which 
clearly indicates that the calculated profile is independ- 
eutec: theslencth of the subinterval [t, _y> thd. The 
smallest time interval entails the most accurate solutions 
and the most expensive run, speaking of computer money. 

It was found that convergence of the numerical solution 
was reasonably fast, accurate and inexpensive, when using 
Peoimeesubinterval Oral min. Theretores in all subsequent 
optimization runs with the DCM algorithm, a 1 min discret- 
ization interval was employed. 

TheSinitial guess of the temperature history T(t) 
ETequlnedw@ to start. themgDCM alcorithm was chosen to be an 
TSG feria leprotile, rer convenience Teasons. me wcood 
initial guess is a matter of good knowledge of the 


physical behavior of the polymerization system. A 
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pertubation of 3°C in the initially guessed isothermal 
profile could cause some serious convergence problems in 
thes DCM alcori thm. 

Ihemecholeesoratne parameters in equation (1iVs12) ais 
challenging, because a very large value of a produces an 
oscillatory or diverging numerical solution, while a very 
small value of a may increase considerably the number of 
iterations needed for convergence. The numerical value of 
Opis directly related! togthe order of magnitude on Gee 
and aGey) termes Tt takes values in the range of #303°K 


6 9 


to 373°K, and g~ varies from 107° to 107 


5 


invehveescudyy. 


a value of a = 1.x10 ~ has been proven to yield a satis- 


factory convergence rate. A value of a equal to o = 5.x107? 


Teads to oscillations, while a value of a = asl 
increases the number of iterations needed for convergence 
from 18 to 41 iterations. The effect of different values 
of a on the optimal temperature profile is shown in 
Pigore 1V°4." it should be noted that’ a chance inio 
changes only the number of iterations needed for conver- 
Sencerand tot the optamal@protile, 

Picecellabulitvyeet scue =DUlisalCOritnneand = tsomconve r= 
Bence schavacleristicspis again®*chnecked ane Figure 1.5. 
Wihebewt NeseGLrect. Or thie snumber OreLterati ones Onethesoptama | 
temperature is investigated. In all cases, an isothermal 
beuiperaluleaproll leva, (UC el SPusecce acm Unem ida Clmpies cm 


As"the number of iterations increases, the profile gets 


steeper and finally converges to the optimal profile. 
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TENE SI Shooting Method Algorithm 


A slightly different formulation of the minimum time 
problem is used in this approach. The objective functional 


Gut sederinedmacmtollows. 
ic 


f£ 
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is 
O 


Final conditions are imposed on Y(t,), as each of the 
fina estate: variables Yo(tg), Y3(te), Y,(t¢) must be 
equa arOnONe mmc lueChismcase G (Y(t,), te) ="0> That means 
Ea ienvem lina scondtl1 one 1V.45 ) for P(te) must change to 
Bie weOld isc Gna (Gly 41 6, ie 

Pi (te) = 0. 

Po(te), Pa(te), P, (te) are free COes?) 


The necessary condition on the Hamiltonian for optimality 


remains unchanged, that is: 


(el Veewl9 ) 
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Dip leis constrained, relationship (21V.19) breaks on the 
boundary and is replaced by the following relationship: 


oan ig, T t) 2 HCY, pe T,t) (IV .20) 


for sail admissible temperatures fT and) all ta [tort,l. 
ASmile@ule DLEVLOUS SOCLELONNOu LOIS chapceim. we nave to 
solve a two-point boundary value problem, which consists 


of eight differential equations (four state equations, 


(vf. VT} 


ajit (ig 7Reo es ee 


$c) Stn (,4),4 ‘ty tige yng mT» 
aneem SANT «0 = (44 tas D.-iige eidi ol . 906 os i 
ey aged samp Yh Sat! ct tl fate bend ota , 
EE 89 aaa Dna: . 


6 = byte: 
et. Vi) 4373 S18 ie feat tpi 
uo ilania'o tt SS sit ned sit ee 

ak 7 bsamuebsiay 


(Ph. 72) ong 
wiih cs enoend (h 12) oa me oF 


gidestijeis: ghfwnllog sia 
Cot. 71) (9,0 wey 
Ages.f) 24 Mean F ” 
as avad. 3% targsdo ‘atts te 


61 


four costate equations) with four initial values 
CuO ae) and four final values (Pj (te) =0.; Yo (te) = 
Y(t) = Y,(t_) = 1.). The Hamiltonian must satisfy 
equation (Volo) forwequation =@ly 20) along the optimal 
trajectory. A shooting method is used this time to solve 
Ene sIRBV problem. 

We shall now outline the steps involved in the 
shooting method algorithm. 
Step 1: Guess an initial value P'(o) for the costate 
vector. 
Steppe eintegrate worwards from E, tO te together state 
andecostate equations (BV.1) and (1V.4) with initial 
condition (1V.2) and the assumed initial value of ey). 
Simultaneously determine the temperature hiEstoryeby 
appryite Trelacions (fv 19) "or" (1Ve20) walone™ the state and 
costate” trajectories, 
Peery 8At final Time (fixed im vadvance)) compare 
Yo(t_), Y3(te), Y, (te) and Pi(te,) with their respective 
desired values. If the stopping criterion described in 
G@quatiom (LV.14) is satisfied, them save the values of 
OUeiMame COMNtLCOM LE theystoppingecniterionsisenct 
Salusi ied @ SO to step 4, 
Steps oe COMpuLeinew initialavaluesstor wthescostate, var- 
iables, and return to step 2. 
ine of aleorithm is summarized In hagunely 26% Ns Ane the 


previous algorithm, the final time is fixed. The minimum 


time is found by solving the TPBV problem for different 
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final times and selecting the smallest one for which conver- 
gence can be obtained. 

Step 4 presents the greatest numerical difficulties 
because, during this Step, s.OuUGetinalevalues (Pi (te), 
Yo(te), Y3(te), Y,(t¢)) must be computed using the shooting 
method. In this step the roots of the vector function E 
relatively to P(o) are determined. F is defined in 


Cqugeloneclye? 1) 


Fy Yo (cae 

pee =|" 24) Sea (IV.21) 
F, a (t,)-1. 
Fy PL (t-)-0. 


A Newton-Raphson method has been used to compute the new 
initial values of the costate vector. 


BO eee sO) oy anne Dene (INE 22) 


Wieleeonicecaurclaxation factors (Og<@¢ucel)., and Diissthe 
(4x4) Jacobian matrix of F relative to EGO ea ihe 


coetficients of D are defined as: 


oF. 
= = i= rve2 
: jalp2pe.4 


ESCs snothantexplicicetunctLlonvor PCojeetheretore the 
Jacobian D must be computed by a finite difference method. 
inaseinvolves four®additional integrations of the state 


and costatevequations and usScarrred*out as follows: 
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(i) Pertube consecutively each component of the 
costate vector P(o), one at the time, while keeping constant 
the numerical values of the other components. Obtain four 
additional initial costate vectors ye Paco PaCcoe 
EaCOn 

BO = ea a Bega 


~ Nee ee 
Ane Ey SACHA Os Daemons sean ys (IV.24) 


(Gieintecratestorwardsmtnemstatesand seOstaiend i cer— 
ential equations jusing initial condition (1V.2) for .the 
Stace vector andythe initial costate vector P+ (0) as 
calculated above. Optimize simultaneously the Hamiltonian 
along the state and costate trajectories. Four new 

3 4 5 


EUNCEVON WWVeCLOrs me, pea yt Fee remthus Obtainedatnat 


correspond to the costate vectors oa, P>(o), eel, Piao jis 


(iii) Calculate the partial derivatives in the 


Jaco manuiaemix AUS ing! a shinvte difference approximation: 


2 ner alebes! gee iol (1V.25) 
Le eee) te CO) ena, 

The pertubation 6 P. was taken as 10% of the corresponding 

Value of the initial costate variable P, (o) during the 
first ten or twenty iterations. As the solution approached 
the optimal trajectory, s P. was reduced to 5% of the 
corresponding value of P', (0). 


The inversion of the Jacobian obtained numerically 


was performed with a very accurate algorithm proposed by 
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Moyer (1978). To avoid singularity of the Jacobian, the state 
variables have to be dimensionless and varying from 0 to l. 
The relaxation factor p was taken as 0.1 during the 
first ten Omestwenty iterations, mcnens it wasechanced) to.0.01 
as the solution approached the optimal trajectory. 
The same stopping criterion as the one used in the DCM 
algorithm was employed to terminate the algorithm, that is: 


RY Ct ee yee t= 2 354 (IV. 26) 


tomadd.LroneesO #LNe sabDoVver condition. Pi (ts) Shouldy pe. Zero 
according to equation (GLV218)i0) (Indeed vit was tound that 

P, (te) was approaching zero as the solution came nearer to 
the optimal temperature profile. We choose to impose final 
@ondLtzons on ¥(t¢) instead ofsminimizine the cosm runctronal 
weormequacion .(1V53) as we did) with DEM aleord thm, because 

Ise waseeLOUNnG that it was Gasiersto Ssatushy, these final  scate 
conditions than to satisfy the Vay incest iia ecOndd tlons som 
P(t¢) imposed in equation (IV.5). 

Tnesinteccation ol the didberentia le equations (lvl) 
and (1V.4) was performed by a Runge-Kunta routine, with an 
AGsUscaDLewineeeratsony SLOPE sazen elem ad Glodc Copies ize 
was five seconds. The differential system did not show any 
stiffness, due mainly to the use of dimensionless variables. 

The optimal temperature profile was continuously 
calculated by calling the minimization routine at each 
integration point. Upper and lower bounds were imposed 


on the polymerization temperature (313°K and 373°K). The 
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temperature was computed by minimizing the Hamiltonian H. 
A half interval method was used to solve equation (IV.19) 
together with the analytical expression (IV.16) for 
(cot Peel eenoeso LUEtONe: Om equatlons( LV. L9) could ibe 
foundsin sthesinterval, (3133/3 )e8the optimal) temperature 


was calculated from equation (IV.20), 
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It was verified that the Hamiltonian was an unimodal 
functron ine chesinterval [3433/3510 meaning that the 
calculated optimal temperature policy was unique in that 
temperature range. 

Convergence of the SM algorithm is highly dependent 
CamLneomintttalectiess of themcostate vector 1 Xo). Since 
the costate variables have no physical meaning, it is rather 
difficult to know even their order of magnitude. To obtain 
reasonable initial values of the costate vector, we used 
the following method: 

Gia rer cube Yo (te), Y3(te), Y¥,(t ¢) around their 

desired values. Set Pj (te) to Zero. 

(ii) Compute Py (t,)» Pa(te), P, (te) from equation 
GiV2o)o) -ASsumera reasonable value for ete), 

(iii) Integrate backwards state and costate differential 
equations (1V.1) and CiVec atc om te to t)- Compute 
simultaneously the polymerization temperature by using 


equations (IV.19) haat (ON PAO 
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Guess four initial values for 


costate variables 


Sa 


Integrate forward simultaneously 


State and costate equations, while 


computing T(t) by minimizing H 


At a 
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Gaioy 788 t, = 0, check the initial conditions on Y. 
Va Y, (0) Paolo sem CNesCOLTesDOnd Np ek( Oymeaamtnel ing tial 
SUCS SE Veo tops mOlatne.oM aleoori thm Lt Y5 (0) is not 


Zero, restart: from (i) with®afatbbterent final #eime-. 


IV.4 Comparison of DCM and SM Algorithms 
The two algorithms were tested by solving the same 
minimum time problem of obtaining a polymer with the 


desired MWD and final conversion as given in Table IV.1. 


TABLE SIVA L 


Desired Propenueies 


fhestesults obtained by the twepaleorithms arespresented 


te Our Ve. 


The difference between the final values obtained by the 
two algorithms can explain the small discrepancy between 
the two optimal temperature profiles plotted in Picure (LV. /)s 
Furthermore, as it has been shown in Chapter II, Blaley “eyeyl 
merization system is not very sensitive to temperature 
Perturbations, at Lhe end of the polymerization. Thus 
the rise in temperature as computed by the SM algorithm 
Pofemnon result in a Ssieniticentechanzesolecoc MuD: 

Let us now examine how both algorithms satisfy the 
necessary conditions for optimality. Both algorithms 


S2creny tie differential equations (IV.1) and (IV.4) with 
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Table IV.2 


Comparison of the Results of the two Algorithms 


DCM Algorithm SM Algorithm 


Yo (te) = X/X, Oe 38 O79b2 


Conversionyx+ OScoa5 OF 6D7 
Y3 (tg) = Wotan 0.998 0.996 


ieee value (1V.2), In the DCM alcorithm: ree) is forced 
pomeactetyathe boundary condition (lV. >) sheuerad@enr 

of the Hamiltonian is continuously approaching to zero while 
the Hamiltonian remains constant along the state-costate 
trajectories, which is a strong indication that the computed 
temperature policy is close to the optimum temperature 
profile. In SM algorithn, Yo(t,), Y¥3 (tg) and ¥, (te) satisfy 
the conditions imposed on them within 5%. The gradient of 
the Hamiltonian is always zero, and the Hamiltonian is fairly 


constant, which indicates that the temperature policy computed 
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by this algorithm is also very close to the optimal tem- 
perature policy. 

For both algorithms, convergence is highly dependent 
On initial guesses. In the DCM alsorithm, an initial tem- 
perature history has to be assumed over the polymerization 


interval [t te] before starting the alcorithm,  lnethe 


oe 
SM algorithm only the values of the costate vector (Omani st 
be known. However it was more difficult to guess the values 
for the four initial costate variables than to guess the whole 
Eemperature (profile. «This is due to the fact that tempera- 
ture has a physical meaning and therefore good initial esti- 
mates can be easily found, especially if some knowledge about 
EnesbDehnavilour of ithe system is available. ‘Costate variables 
have no physical meaning, therefore it may be very difficult 
Com indesome reasonable values tor 2 (5), 

Twenty to forty iterations are usually needed to obtain 
Converzence with both algorithms, depending ony thesaccuracy 
of the initial guesses, and on the relaxation factors a and 
op. For the DCM algorithm, one iteration involves the inte- 
gration of two systems of four differential equations. 
For the SM algorithm, one iteration involves the integration 
of five times eight differential equations, and the solu- 
tion of the equation (3H/3T) = 0) at all integration points 
Cisualivel?, 000Spoints). Konecompariconsieacons po tia 
algorithms were implemented on a very slow computer (HP 1000). 


Ten to twelve hours were needed to get convergence with the 


DeMealeorithm, against three sor Tour daysoror (chessiealzo- 
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rithm. The same test performed on the Amdahl 470/V compu- 
ter at the University of Alberta gives a CPU time for the 
DCM algorithm of approximately one humdred seconds. 

The fact that both algorithms predict the same optimal 
temperature profile using two different numerical methods, 
that the gradient of the Hamiltonian relative to the 
temperature is close to zero, and that the Hamiltonian is 
almost constant during the time of the reaction, strongly 
indicates that the solution obtained by these algorithms 
is close to the true optimal profile. However, it is 
difficult to show in a mathematically rigorous manner that 
the computed solution is the true optimal. 

The DCM algorithm has two main advantages over the SM 
algorithm: (a)” it nas Laster convergence rate, Coe iter s 
more efficient. Therefore, for computational reasons, all 
optimal temperature profiles subsequently presented are 


computed by the discrete control method algorithm. 
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CHAPTER V 


OPTIMAL 2 EMPERATURES POUT CLES 


V.1 Reduction of the Batch Polymerization Time. 
The isothermal polymerization of MMA at 50°C requires 


Sereacce One 1MemOLe42 min Porattain: a conversion of 90). 
The M, and Me of the polymer product at the end of the iso- 
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Table V.1 
Final MWD for the Isothermal Polymerization 


OL MMAMat UG 


M M PD u 
n W O 
6 7 -4 5 
: CF IG L spa0 0.406 x 10 Cee WGe 4 5 exeLO BS399x= 10 


Using “the DCM *alzorirehm *to calculate the optina M*temperatire 


policy, we were able to obtain a polymer product with the 
BpeCerticaerors Of tap le Vs only rinte2 > Tiina Deter Sreawe 
gee 25.9, decrease in the “total polymerszarromet ime *com- 
paced tomtie =1sOthermalerrin at750° Greeihe calcu Pared) opirme! 
temperature profile is shown in Figure V.1. 

The isothermal “polymerization ob@MMA*at 870° Créquines 
avreaction time of Sl "min *to“reach *a-conversion om 9072 @ihe 


corresponding MWD obtained at the end of this run is presented 


in Labbe Vr2. 
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Tapireny 2 


Final MWD for the Isothermal Polymerization 


Ont MMAsiaita7.0 FC 
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Using the DCM algorithm to solve the minimum time problem 


’ 


we obtain a polymer product with the specifications of 
Table V.2 in only 67 min. The batch time is therefore 
reduced by 19%. The optimal temperature profile is shown 
iInerreure V2. 

path optimal profiles (Fieures V.i sand V.2)) have the 
Same shape. [Up to about /0% conversion, the temperature 
remains almost constant, then it sharply decreases and 
finally, at about 90% conversion, the temperature increases 
to a final value slightly above the corresponding isothermal 
temperature. Under isothermal polymerization conditions, the 
factesor. polymerization would sharply ancrease as a tesult 
Stmciceceimertect (Chapter Ll). Tniseing urn would, causera 
large increase in wy due to the decrease of the termination 
Bacescoustaie.. Li should be notedsthnatsthe cel eliect is 
more important at low polymerization temperatures than at 
higher ones. Thus if we wish to increase the polydispersity 
ey) of the final product, a decrease in the temperature 
profile at the onset of the gel effect phenomenon would cause 


a dramatic increase in 7m) and a subsequent increase in the 
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polydispersity of the product. On the other hand, if we 
wish to obtain a product of low polydispersity, the tem- 
perature profile should increase to suppress an increase 
in Uo due to the onset of the gel effect. In Figures IV.3 
to IV.6, the monomer conversion and the rate of polymeri- 
gation aresplotied as"a function om reaction time.) It can 
be seen that in this case, the optimal temperature policy 


has a tendency to flatten the rate of polymerization profile. 


V22. Quality Improvement of the Polymer (Product 


The minimum time problem has been solved to produce 
polymers with different molecular weight distributions. 
The corresponding physical properties of a product obtained 
under these optimal temperature policies cannot be gotten 
by any isothermal run. 

By fixing the number average molecular weight Mo to 
eave memequs): to that obtained atsthe end jor gthe 50°C pf iso- 
thermal run (™, = 991,000), and byivarying the weight average 
molecular weight Mo different desired polydispersities have 
been obtained using the DCM algorithm. The desired MWDs 
and the corresponding computed optimal values are presented 
in Table V.3. Figure V.7 shows the computed optimal profiles 
for products with different polydispersities. The correspon- 
ding profiles for the rate of polymerization, M,, M, and 


polydispersity are displayed as a function of time in 
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Valve equalgrogtnat obtainedsatethebend of thet 502 Cn ea soe 
thermal run (M = 4.06 x 10) and by varying M,» another 
SCE OLeOpLimaiyprofilesshas been obtained. The different 
desired MWDs and the corresponding computed optimal values 
aeeeeeven singlablesVvedues The corresponding optimal tempera- 
ture profiles are presented in HieuregV. lives Thesratesot 
polymerization, Mo Mh and the polydispersity associated 

i PEN@Eheceeoptimaimnrotales areashounadn Figs ese Vel? eie13 
andeVel 4. 

When the minimum time problem is solved to obtain 
narrow molecular weight distributions (polydispersityee Z2ie 
3.0) the optimal temperature profile shows a pike-shaped 
increase in temperature around 70 - 90% CORVELSAOR Shox 
Ehescase of broader MWDs (PD, = belies DembeO).. the computed 
optimal temperature profiles show a pike-shaped decrease in 
temperature around the same conversion range. When a polymer 
PEOdUCte Ofehighepolydispersitys (PDs= SU 0S 620) Sa sPdecareu. 
the rate of polymerization shows SLLONS Variations and even 
a bimodal form. As shown in Figures V.10 and ViL4, the 
polydispersity always goes through a maximum at around 
60-60% conversion. “This is due to the acceleration in the 
rate of change of np, duegiic the teclferte.t tees ceerpangce, 
conversion, the polydispersity of the polymer drops dramati- 
cally as a result of the decrease in (Ky = esse hia which 
involves a limiting Wo and ay Jimitine= Conversions 

By Comparing the” two Sets of Optimal proriles shown 


in Figures V.7 and V.11, it can be seen that shorter minimum 


batch times are obtained in the second set (where We is fixed 
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and M fs*varyine) . From :Tables4V.oeand V.4. ®itvcan*be 
seen that in fact the shortestminimum time is always obtained 
for the smallest desired Mo: 

From the™above*results, “it®is expected ™that @any sensible 
molecular weight distribution can be obtained by computing 
the optimal temperature policy using the discrete control 
algorithm. 

As shown in Chapter II, the system responses are sen- 
sitive to variations in kinetic parameters. Therefore it 
may be anticipated that perturbations in the kinetic para- 
meters will result in some changes in the minimum batch time 
Sami lea cneCmOpecimalatemperatunre: proti ver foerind soled © pany 
perturbation in the kinetic parameters affects the shape 
of the calculated optimal temperature profile, the minimum 
time problem was solved several times, each time perturbating 
one kinetic parameter from its nominal value. Results of 


Eiiomanalysis are, Summarized in lable Vn) cand in) hiciresVeloe 


Tabi eave 
BPerect ol barameter Perturbations onmsuneltlin imum 


Lime eh eobeen Problem 


Computed 
Minimum 67mn |67mn |67mn |65mn | 65mn |65mn |69mn |67mn 


Time 


Parameter 
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The kinetic parameters Ay; Ay; Bi, By have a significant 
effect on the optimal temperature policy and on the com- 
puted minimum time. This is understandable since the con- 
version and the second moment are exponentially dependent 
on Ay: Ap: Bi» By Parameters. Again, it can be seen “that 


much care is recommended when estimating these parameters. 
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Figure V.8 
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CHAPTER VI 
OPTIMAL GATALYST “FEEDSRATESPOLEGLES 


Ups toe ehis= point =the only controlgvariable considered 
was the polymerization temperature. In this chapter, 
the minimum time problem is solved by optimizing both tem- 


perature and catalyst feed rate. 


Vi.l Problem Formulation 


Let ies be the catalyst feed rate. The polymerization 
OreMMATing a batch reactor 1s described by thesser or 


differential equations II1.11 derived in, Chapter II1. 


T.. #) CA desis) 


The Hamiltonian of the system can be expressed as follows 


(eq mL t. 16) 
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F. appears linearly in the Hamiltonian, thus Lhe necessary 

G 
condition for optimality (dH/dF, = 0) does not allow us to 
calculate the optimal policy for F,. Therefore the Pontryagin 


Minimum Principle has to be applied to calculate the optimal 


feed rate take 
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for all admissible a ang aller [t); ty]. 


Design reasons require that aoe constrained between a 


lower and upper Hig value, 
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Ce ce max (V1.4) 
Taking into account the constraints on ae equation VI.3 can 
providesthe solution \to,.thefoptimal catalyst feed rate policy. 


This policy will be of the on-off type given as: 


TEP. a> 60 le Ve 
e 
lige Py < FO ne = EA ne CESS) 
ilas Py = oO ; ae is undetermined and can take any 
value between o and he Ls Suchifan loptinal tpolticwerstcalled 


bang-bang control. An attempt to solve the minimum time 
problem by applying this bang-bang policy to our numerical 
model did not produce any meaningful results due mainly to 
convergence reasons. This may be due to the fact that EOu 
eertOcesmaul on too Large Ee oe value, a SOREN Ses to the 
minimum time problem might not exist. In order to find the 
Optimal policy, [I.] and De aes should be optimized simul- 


taneously with the computation of the optimal catalyst feed 


Lace policy. An alternative is to consider 7 lie etnekinitiacor 
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Conecentration,»; asethe control variables instead of iD et 

[I] can be manipulated by controlling the catalyst or inhibi- 
tor addition rate ntoythe reactor. Once the optimal profile 
for catalyst concentration is known, the catalyst and inhi- 
bitor feed rates can be calculated from the mass balances 


on catalyst and inhibitor species: 


Seas = CP] Pee (Nit “e (VI. 6) 
d[N 
GLE} =~ ky II[N] + a CE) 


wheres [N] is" the"ianhibitor= concentration in®the* reactor , 

Fy moe the anhibrter feed rate, and ky Usmunemrate constane 

freconcumptLonmor Iii tiatorevadicalssby thesginhib itor: 

Tt) is assumed. that’ there’ ais’ no? change in the volume of) che 

feactine mixture due to “catalystvand inhibitor saddicion- 
Thespolymerization can bevdescribed sbyatureesstate 

variables (dimensionless conversion, dimensionless yee” 


moment, dimensionless second moment) and two control varia- 


bles (polymerization temperature, initiator concentration) - 
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The Hamiltonian H' for the above system can be expressed 
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wer ces Loe Optimally initiatorgrolacy 


In this case; we consider ithe initiator concentration 
uli the=reactox as{ithe onlyscontrow tvariable. ® Ihe 
Eemperatune © @ss| arbitrarily faxed jtoway constant walue. 
Thesalgorithm used toysolve thus problemiis™similare to the 
SleCuienmensedminecChapter V tomcalcularesthe yomtimaiy tem 
Perature procive. ihe initiator concentrabvionmls. updated in 
each interval [ty _1> ty] DY Applying sa sist mOLdermme madden 
method: 
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magnitude as the term [I]. A value of a = 2 was proven to 

yield a good convergence. The gradient of the Hamiltonian 

relative to the initiator concentration can be calculated 

from the following analytical expression: 

oH -¥ 32 | Negus pi + praca. pi (VI 15) 

d[ 1] Zi) Xa Hod 2 ‘ 

The numerical solution starts by assuming a constant profile 

FOL [I], throughout the polymerization time from t, to te. 

Again the “discretization time “interval was chosen as 1 min. 
The minimum time problem has been solved for two desired 

molecu am wetehi=tdastrmbuUtions, “siven iim Pabple Vil” The 

corresponding optimal profiles are shown in Figures VI.1 

and VI.2. For the same desired MWD, the shape of the 

Opiamal catalyst concentration prorile is “Similar to the 

shape of the optimal temperature profile. This is under- 

Stancepley ‘since’an “increase in “initiator concentration and 

an increase in polymerization temperature have a similar 

erfect onthe state variables.“ as™ shown on Figures I1.12 

and II1.13, where temperature and [Io] sensitivity coefficients 


were presented. 


Vinge thesOptimal Temperature andj initiator rolicies 


In this section, both the polymerization temperature T 
andmthencatal yst concentration) [Ely vavesconsidercdges 
control variables. The DCM algorithm is used to solve the 


minimum time problem. The two control variables are updated 
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using the following equation: 


3H 
UNSR mew amc old oi om le, 
(VI.17) 
Q oH 
CONG an) EMIS Oe Hay 


Paemnumericalevaluessol O47 and a, were Chosenmacmuncc 10° 


and 2. respectively. For these values, reasonably fast 
convergence was obtained. The initial guesses swepe (Gm) shave! 
[Il(t)were respectively an isothermal profile (50°C) and 


ly. 


The optimal temperature and initiator policies that 


a constant initiator concentration (0.02 mole.2 


minimize the batch time were calculated for two desired 
MWDs described in Table VI.2. The optimal initiator concen- 
tration policies are shown in Figures V1.1 and Vaz and 
compared with the optimal initiator concentration policies 
computed by considering only [I] as a control variable. 

The optimal temperature policies are shown Ta Pieure sa vd ie 
and V1.4 and compared with the optimal temperature 

policies computed by considering only T as a control variable. 
The shape of the optimal temperature profiles and of the op- 
timal initiator profiles are alike for the same desired MWD, 
but the simultaneous optimization of [I] and T has a smooth- 


ing effect on the profiles. 
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VI.4 Controllability of the System 


The system of differential equations which describe 
the polymerization of methyl methacrylate is non-linear, 
therefore it is difficult to show or to demonstrate if the 
system is "controllable" in the technical sense ,94i Je, if 1c 
is possible to transfer an arbitrary initial state to an 
arbitrary final state in a finite time. In this work, con- 
vergence is always obtained because the final desired out- 
put values are allowed to lie anywhere within a +5% band, 
Which is satistactorys from jaspract ical: engineering point 
of view. If a smaller tolerance were imposed on the final 
desired output), it mayshavexbeen moretdiffaculteand, in 
some cases, impossible to obtain a true optimal control 
policy. This observation is supported by the following 
results: 

For two computer runs, the desired final conversion, 
Oth and 2nd movements of the molecular weight distribution 


were specified as: 
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For the ees run, a tolerance of +5% was imposed on 

the final desired outputs. As shown in Table V1.3, con- 
vergence was obtained in a reaction time of 290 minutes. 
For the second run, a tolerance of +1% was imposed on 

the desired final output. Even by increasing the reaction 


time to 335 minutes, the difference between desired output 
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and the actual output was reduced only to 3%. 


Table.V1.35 


Effect of Tolerance Limits on the Convergence 


As it was shown in Chapter II, at the end of the reaction, 


Erroronen 
WwW 


the conversion and second moments are not very sensitive 

to temperature and initiator concentration (the two control 
variables ye ethisemeans that iteis di cercu kt COsgerearic 
desired final value for conversion and second moment once 
we have passed a critical point after which the sensitivity 
coefficients for conversion and second moment get opposite 


signs. 
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CHAPTER VII 
CONCLUSION 

A mathematical model has been developed for the free 
radical bulk polymerization of methyl methacrylate initiated 
by AIBN. The dynamics of the polymerization system are 
described in terms of monomer conversion, zeroth and second 
moments of the molecular weight distribution. The kinetic 
parameters of the model have been estimated by fitting the 
experimental data of Balke (1972). The sensitivity 
coefficients of the system outputs relative to the 
kinetic parameters have been calculated. It is shown that 
the second moment of the molecular weight distribution 
and the monomer conversion are very sensitive to the poly- 
merization temperature, and to the kinetic parameters Ay» 
By. A>» B., describing the autoacceleration of the reaction 
due to the gel effect. Therefore much care is recommended 
in the estimation of the above parameters. 

Our main objective of applying the optimal control 
theory to a polymerization system has been accomplished. 
The minimum batch time problem for the bulk polymerization 
of MMA has been solved by considering the polymerization 
temperature and/or initiator concentration as control 
variables. Using the calculated optimal polteLres we were 
able to reduce the batch time of the corresponding isothermal 
batch time by 25% while keeping the same polymer quality. 


Furthermore, we show that a desired MWD can be obtained 
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by applying the computed optimal temperature and initiator 


concentration policies to the polymerization reactor. Two 


algorithms (DCM and SM) are developed in this work to 


solve the two-point boundary value problem generated by 


the application of the Pontryagin Minimum Principle tosthe 


mathematical model of the polymerization of MMA. Although 


both algorithms give comparable results in terms of minimum 


batch time and optimal temperature profiles, the discrete 


Control faleori thm has proven) to be the: bastecteandmmoce 


reliable. 


Asva result lor vthis investigatiom, the etollowine suture 


projects are suggested: 
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C2) 


(iii) 


Experimental wverification, of sthe optimal policies. 


Optimal policies and particularly the optimal 
temperature policy can be implemented easily on a 
DiLOUypoOlymeni zation Teacvon » sUCh 45D. orp lant 
system has been built, and is at this time ready for 
implementing time optimal control algorithms. 
Extension ob this "work ito others chemica lr ecacuionse 
The algorithms developed in this thesis are very 
general) and can, be) used for the optimization Osa 
large number of chemical reactors. Reactor 
productivity couldetherer ores Desecrea tl amo roved: 
Reduction of energy consumption. 

Other objective functionals than minimum batch 


time should be considered. The expenditure of 
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energy to control the reactor temperature should 
be taken into account. Still the same algorithms 


could be used by deriving an other expression for 


the Hamiltonian. 
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APPENDIX A 
DERIVATION OF KINETIC EQUATIONS 


The mechanism of the bulk free radical polymerization of MMA 


initiated by a free radical catalyst is described as follows: 


Catalyst decomposition 
ik cRewtatust: eporense 


Initiation 


k * 


ee ae 
Be tat . Ry 


Propagation 

* 

Feo. M Be eee R for xen lt, &) 
Xx X+] 


Termination by disproportionation 
> Eee 
R, + R, —__ft _ "x *./y for (x, y).e [1,.¢) X [1,.%) 


Where I represents the initiator molecules, Ro the initiator radicals and 
R. the activated initiator radicals. M denotes the MMA monomer. R is 
the live radicals of chain length x. ane js the dead polymer of chain 
lengthax, «of siSs the initiator efficiency factor. Ky kK. ky k, are the 


rate constants for initiator dissociation, initiator, propagation and 


termination reactions. It is assumed that rate constants are independent 
of chain lengths. 
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According to the above kinetic scheme, the mass balance equations 


for the different species in the reactor are derived. 


Catalyst 


The mass balance for the catalyst can be expressed as: 


at eps 
toa Ky eo, Aue 


Where [I] is the catalyst concentration in the reactor, and t is the time. 
If catalyst is added during the reaction, it is assumed that the quantity 
of catalyst added will be small, so there is no change in the volume of 
the reacting mixture. In this case, the mass balance equation for catalyst 


is written: 


sl J nate sty 


Where i is the catalyst feed rate, and V the volume of the reactor. 


Monomer 
Monomer is consumed through propagation and initiation reactions. Gener- 
ation of monomer results from termination reaction between a free radical 


of chain length one, and any other radical. 


co jee) 
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d [M 
dt : : x=] x=] 


* 4 ‘ F 
where [M] is the monomer concentration, [RJ is the concentration of live 


radicals of chain length x, [R is the concentration of catalyst activated 


radicals. 
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Activated radicals of catalyst 


* 
Gol = OF k 


[ L * 
ae ete altel) (4) 


Activated radicals of chain length one 
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1 =e * a ) 
m K. [M] [Ro] ky, [M] [Ry] k, [Rj] [RJ (A.5) 
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Activated radicals of chain length greater than one 


co 


alee Ri (ML Sebk TR elie k [R.] [RI 
= pen x Diex= 1 TOieR Yor itAse ) 
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vi mes [Rel y ) [pl (A.7) 
dt co; y 


Where Bis is the concentration of dead polymer of chain length x. 


The moments of the MWD are defined as follows: 


yn moment of live radicals Gistripuuion. 
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The rate of change for the moments of the MWD can be obtained by 


differentiating equations (A.8) and (A.9). 


Oe) owe Sy med [RJ cert 
dt dt ; 
xX = ] 
d 
d ies SS F d BRS | 
dt ‘s dt (A.11) 
Kt= or? 


Substituting equations (A.5) and (A.6) into equation (A.10), a 


differential equation for the yth amount of live radicals can be defined: 


[oe) 
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Using the definition of the nat moment of the live radical distribution 


(eg (A.8)), equation (A.12) can be rewritten as: 
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In order to simplify the above equations, the quasi steady state (QSS) 
approximation for the live radicals is employed, which means, that the 
derivatives in equations (A.4), (A.5), (A.6) can be set to zero. Similarly 
the derivatives in equations (A.14), (A.15) and (A.16) can also be set to 


zero. Thus equations (A.4) and (A.14) can be rewritten as: 


* 
Peale oka (MI (Re 1 (A.17) 
* 
on op Kay tM) CBG! (A.18) 
eres 


Combining equations, (A.17) and (A.18), an algebraic expression for hy can 


be obtained. 


Using the QSS approximation, equation (A.15) can be recast into 


equation (A.20) 
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Since the term (2f ky [I] ) is negligible compared to (k, [M] A, ), the 


following expression for A, is obtained: 


' (A.22) 


Justification for these approximations is given at the end of this 
appendix. 


Similarly, equation (A.16) can be rewritten as: 


hy = 2f Ky Ble 2 ky [M] Le ky [M] Ns 
k 


(Avo) 
t 


Since O09) is negligible compared to (a4) and (2f kg [1] ) is negligible 
compared to (k_ [M] ho ), a simple expression for do may be derived from 


Pp 
equation (A.23): 


=1 2 k [M] hy 


Ke do 


(A.24) 


Substituting equation (A.7) into equation (A.11), differential equations 
describing the rate of change of the zeroth, first and second moments of 


the dead polymer distribution may be derived: 


. 2 

dig = ky Ag (A.25) 
dt 

dup = eo 25 (A.26) 
dt 

dug = Rtg fe (A.27) 
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Substituting equations (A.19), (A.22) and (A.24) into equations (A.25), 
(A.26), (A.27), we obtain: 


du, 

ay gp ak LD (A.28) 

vi 

ele. far ee BS, 
dee amy, My a (A.29) 

t 

d_, ee 

os Aoi See (A. 30) 
dt 2 — [M] 


The monomer conversion X is defined by the following equation: 
Neel rane Ie 
fo) 


Where [Mo] is the concentration of monomer at time t=o. By differentiating 


equation (A.31), the rate of change of conversion is obtained: 


dX ] d[M] (A. 32) 
dt 
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dt [Mo] 
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Since (2f ky Gh Se k, [Rj] ho) is negligible compared to (ky [M] ods 
equation (A.33) can be simplified. 


; (A.34) 
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Substituting equations (A.19) and (A.31) into equations (A.34), the 


following differential equation for the conversion is derived: 


We Bets, WU eS te) (A.35) 
ae 
k 


By noting that the sum of the first moments of the dead and live polymer 


distributions is equal to the monomer consumed: 


tee, [M1 - [M] (A.36) 


Neglecting the contribution of a jin equation (A.36), we obtain: 


Hy ath X (A.37) 


Differentiating the last equation, we get: 


d Wo [M 


dt 


} dX 
C tdt 


This fast result can be verified by combining equations (Ae29 Wanda (A, 35). 


Thus, the polymerization of MMA in a batch reactor can be completely 


described by the set of equations (Ae2\c) (AGS )ee(AC28) mand Aes0)e 


Justification of the Approximations in Eouationsa\ nec, se hecsis ann oon 

This section intends to justify the approximations made in deriving 
equations (A.22), (A.24) and (A.34). Table A.1 shows the numerical values 
of the different terms which appear in the kinetic equations. The 


following numerical values for the kinetic parameters have been taken 


after Sacks et al. (1972): 


b i ed 


(ae) 


alana a (Bei) anes # tno 
_ 7 Nee 


_ 


(ig. a} +f | 5 ee a 
e| =. 
7 


: Sonn sw oOFgupe sent st | 
, Toe oe “pins 
| — a 


ato “as 08,8 nh in eran 
Ufpsalipons ed ee vides ee 


pelea 

fa 

ae re 
teeto’ fyi 


oan 


si in 


Ky Sibel os Ya @ x EXP ( - 30800/RT) so} 


ero 1X 10° x EXP ( - 6300/RT) 1. mole! s~? 
_ 8 iba aa 

k, = 7.8 x 10° EXP ( - 2800/RT) 1. mole” s 

f = 0.6 


l= wal x 107° mole/] 


[M] 5.mole/1 


From the numerical results of Table 2.1, it can be seen that es 
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APPENDIX B 


DERIVATION OF THE SENSITIVITY EQUATIONS 


The bulk polymerization of MMA in a batch reactor is described by a 
set of four differential equations (A.2), (A.35), (A.28), (ARB 00 


dZ Fee pe) (Bey 


and a set of four initial conditions: 
(sea (B.2) 


where Z is the output vector and p the parameter vector. 
- if 
Z is Cea Ks Uo» Uy) (B.3) 
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These parameters have been identified by fitting the experimental data of 
Balke (1972) with our own mathematical model. The purpose of this appendix 
is to describe the steps required to calculate the sensitivity coeffi - 
cients that show the effect of pertubations in the parameters on the output 
of the system. 

The sensitivity coefficient for the parameter Ps and the output Z. 


js defined as the first partial derivative Gh i. Wit hares pect tO PS: 
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The sensitivity matrix is the matrix of the sensitivity coefficients: 
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The information provided by the model must be used to determine 


=I Q 
fo] |™N 


However, Since equation (B.1) cannot be integrated analytically, 

the sensitivity analysis cannot be calculated analytically from equation 
(B.6). By using the sensitivity analysis method (Atherton et al. (1975), 
Beck and Arnold (1977)), a differential equation for the sensitivity 
matrix may be derived. To do so, we interchange the order of did renen— 
tiation in the expression below and use equation (BW )to obtain: 
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ate 3 


Seer aloo FOC CVON OF ee Une right hand expression is expanded to give: 
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The matrix of derivatives are defined in the usual manner as in equation 
(B.6). By introducing the sensitivity coefficient o5 defined as in 
equation (B.5), equation (B.8) may be written in the following form: 

j=4 
“ij afi, oy | (B.9) 
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To determine the initial conditions for equation (B.8), we note that at 
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These conditions are summarized as: 


Omit P, is not an initial condition 
(o) = (B.11) 


oj if oe is an initial condition 


Equations (B.9) and (B.11) give a set of differential equations whose 
solution determines the numerical value of the sensitivity matrix as a 
function of time. The number of equations is equal to the number of 
state variables time the number of parameters (i.e. 4 x 11 = 44). It 
can be seen that each sensitivity equation is linear in the “ae and 
has variable coefficients which are determined by the state Z. 

The equation (B.1) describing the polymerization system has been 


derived in appendix A (eg. (A.2), (A.35), (AP 26)0 UR SO }ineeel te Seas 
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Application 


v2 


of sensitivity equations (B.9) and (B.11) to the system 


(B.12) yields the following differential equations for the sensitivity 


coefficients: 
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For initial initiator concentration [I] 
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The numerical integration of the sensitivity equations has been performed 

along with the numerical integration of tie strate equations, using a fourth 

order Runge-Kunta with variable step Size integration routine. 
Alternatively a finite difference method could be used to compute 

the sensitivity coefficients. The finite difference approximation for the 


aa output variable and the parameter pj is defined as: 


02. : i. (P,> Pos Rae) Ps a OPs> a Pry) 
en dk 6P 5 
j 
6p. 


J 


7 _ _ 
oo vorersind fotatatimel 


/ 
; 
4 
a 4 = \ ; 7) 2 e ¢ : 7 
= ih ~ > (9) tis ‘ (o) [ii 2 
- ea 
ie a 


barwtisq otal zl i tmees yRiwhdbenes sty Yo not davgeret isohyamun ‘af 
dite? # pattu partes odnee Gis Be nat twrgeznt feaitar wn 3c thw ¢ note 
ri bugt HOMetadat abek qege Gidatvey Agrw nzewl-spmull 8 we 
chuqees G) sey. GO Lite Setpew goawpysvith estar? 6 peers 
" gity oP nai teahenhagg Ganenet hie SETPE WAT .enatdiy ted one rte ’ 
a6 Dental. ef nay a 


) ae : 


oe 


3) 


Following the recommendations of Beck and Arnold (1977), De has been 


chosen as: 


Spe e=oeO20 ep. 
PS 0.0001 p. 


This method was used to verify the correct integration of the sensitivity 
differential equations. Sensitivity coefficients calculated by the two 
different methods, did not vary significantly. 

In order to show the relative influence of the parameters on the 
output variables, a normalized sensitivity coefficient was used in place of 


the true sensitivity coefficient. It was defined as follows: 


These normalized coefficients are plotted in Figures I1.8 to i lal 
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APPENDIX C 
DERIVATION OF NECESSARY CONDITIONS FOR OPTIMALITY 


This appendix presents a heuristic proof of the necessary conditions 
for optimality derived from the classical calculus of variations. 

A detailed proof of the necessary conditions for optimality and of 
the Minimum Principle can be found in the work of Pontryagin et al. 
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The process to be optimized is described by a set of differential 


equations for the state variables together with initial values: 
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where Y is an n-diniensional state vector and u is an m-dimensional 


control vector. Suppose that we wish to mimimize the functional J, 
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Where G and L are two scalar functions. 
Introduce the n-dimensional costate vector P (t) and rewrite the 
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Where the superscript T indicates transpose of a matrix. 


In order to minimize J, we must make 6J = 0, 
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Expanding the terms in equation (C.4) to the first order, 
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Substituting equations (C.5), (C.6) and (C.7) into equation (Or4 Sand 


using equation (C.8) to calculate tne first order variation “J, we 


obtain: 
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Sances(f-Y-0) =the coeiiicient of sp! vanishes. 
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: pl sy dt 
Integrating by parts the integral ous we obtain 
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(C210) 


ean 


Substituting equation (C.10) into equation (C.9) and rearranging, 


we get, 


eae} ae Tap ¥o mera Ys09 aily (ost?) sont 


: “4 
i 
an ‘we wa “S \ gana! at? eyreq vd gnttnigazal 
r on 7 ; 
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vgs ea ae 


139 


qa 9) 

6J 3 Pp OV(t ey ss 

te 
_ if t aT I ie RNES 

ine fe at [eet sp. an dt (Cui2)} 
a + oY a ie are Pare i 
ay JU cu 

t 

0 


éY and Su are arbitrary variations, therefore necessary conditions 


for 6J to be zero are that the coefficient of 6Y and du are zero: 


aL of : 

yr ier ac 

Ty Senpe Gee (Cale) 
ee re ok 

aaa — = sae! 
du lee (C.14) 


Define the Hamiltonian H of the system as follows: 


Feoevus = GC the mac (c.15) 


Equations (C.13) and (C.14) can be rewritten in terms of the 


Hamiltonian: 
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J 3H 

p= - 

oF oY (Cis) 
3H 

sua (Cage 


Suppose that Y.(t,) ist ect Isom ece——| 


and ¥5(t,) ie SRRAD Gdelee se = fe hy Wear 5 Sen 


then sYi (t ¢) = Gy sie Wt a2 ele 


and oY (t,) ican arbitrary Variation 101g -shate).--0 


The following final conditions for P(t.) can be then derived from 


equation (C.12) by forcing 6J to become zero: 


aG 
P (te) = aye for j= ret 1s ot 22-——n 
Ji te 
P.(te) is free fOr = Iecseet (Cal3) 


It should be noted that in the minimum time problem, H is not an 


explicit function of t, thus (dH/ot) = 0. 


Differentiating H with respect to time: 


dH 3H 9H du oH dP oH d¥ 
ae ee hy ee ee es | eae Gane ha 


dt.) (ots OU ct oP dU ee (C.19) 
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However, 


along the optimal trajectory, the necessary conditions for 
optimality will be: 


dy 
aD? GR Ape an” 


(C.20) 


Therefore; H iS constant along the optimal trajectory, 


(dH/dt) = 0. (C21 
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APPENDIX D 


DEM and SM Aleorithms Computer Programs 


* 


‘amergort tutiqmod aaiaitogiA ME baa Mod 


FING 

c 

c 

C SHOCTING METHOD ALGORITHM 

c 

c 

C MINIMIZE THE HAMILTONIAN BY 

c 

Co GUESSH PMO» 

C DIMENSIONLESS Y1,Y3,Y4 

CMA CEs SEs Nae teniret, Ys Om Ea, yie 

C INTEGRELTE FORWARDS 
EXTERNAL MINIH,RKF,FCT,O 
DIMENSION Y‘&:;,DERY(8',Y 
REAL KD,KT,KTC,KP,K1,MO, 
DIMENSION GUESS‘5,4),Fi‘5§ 
DOUBLE PRECISION FACT( 4) 
DIMENSION NAME(3),10CB'1 
COMMON TEMP .MO,10,HAM, DH 

c 
WRITE( 1,20) 

20 FORMAT(1XK,‘ENTER FIRST T 

zs’ WHERE THE HAMILTONIAN 
READ(1,#;5)LU1,LUz 

c 

c 
WRITE(1,21) 

21 FORMAT!1XK,‘GIVE FINAL TI 
REACO'1,2 )MFIN 
TFIN=FLOATIMFIN! 

WRITE AY, 231 
READ, FVYO C2) YOUS+:, YB 
WRITE!1,24) 

26 FORMAT: 1X,’GIVE FIRST GU 
RELL tS VUGVESS It) vat, 
WRHINE HEUie 2 200 YbiGNa = 25 

22 FORMAT /1%. “YC2VOESTRED=" 
ieee Wee SmRE Dim oy) E1013) 
WRITE DEURE 25 GUESS ( ) 7 J 

25 FORMAT(1X,'° INITIAL VALUE 

c 
(S 
c 
C INITIAL CONDITIONS 
Mo=c.3 
10= o169s 
Wit ert 
Yu20=o 
v131=0 
YI&r=0 
WES P= GUESS ey 1 
yre Wi v2 
Pe T=IGVES SE iTS 
YREV=SGUESS ¢ 1.8% 
C START INTEGRATION FORWARDS 
NOIM=E 
x=0 
b=) 
BDC 6O K=1,MFIN 
BO 762 Wet, Ne 
DA=+5 
HH=+S 
HMx=S5 
RELER=0.1 
ABSER=0 
GALE ROOTHT), Hor T , TOP T ; 0 
TEMP=TOPT 
HAM=HOPT 
DHOM=DHOPT 
EAP EORKEEX  Y,NDIM, FCT, DA 

G2 CONTINUE 
MM=rk-L 
L=K 
CALELOUTPCK,¥ TFUAG, CO 

60 CONTINUE 
WRIREOEUD PIFEKP OV EI, Pet, 

7 FORMAT T// / 1X, “RESULTS (FO 
PIRGh Meo ek, SEV Fy Lh Ks 

€ TEST FIRST GUESS 

Eis hore ys ha 

Fait eshey ta c= 4 

Rie smevr a) =4 
F(y,@reyrcs' 
BREABSPER YS. tht 

TRGPE OG Oo, O25 Go TO 65 
EFPSARS WE Oi); 210 

IF‘FF GT.6.025 co TC 65 
FREABS LE ht 30} 

HRCRA GY OroO2S 7 eGo TO “és 
FFZABS'(Ft1,4)) 

VRUEF ECT. GeozsS! GO TO 65 
WRITE(LU1,66° 

66 FORMAT:1x%, OBJECTIVE FUN 
GO TO 40 

65 CONTINUE 

c 
c 


USE OF NEWTON-RZLPHSON METHOOC FOR THE 2-POINTS BOUNDARY 


DH/DU=O HALF INTERVAL 


TF) 


UTP, H,DERH,ROOTH 
Dia) PFI 4) 
J£C(2 6: ,306 

,4) , DELTEL1 4) 


METHOD 
VALUE PROBLEM 


DOUBLE PRECISION AMAT(16),A1116+,R116),R1/'16),EPS BB‘ 4: CO(4: 


fe LG. 
44) ,1S1ZE(2) 
am,yYD 


HE UNIT NUMBER WHERE 


INTEGRATION’, 


sVeReESUUTS ARE TO SE PRINTED /IUK ANG THE UNIT NUMBERS, 


HAS TO BE PRINTED’? 


MES TNT EGER) <1 


22 FOPMAT( 1X, 'GIVE DESIRED CONVERSION DESIRED MOMENTS’! 


1COUNT=1 

1000 CONTINUE 

COMPUTE ANC TEST FOUR PERTUB 
DEL= 1 
FaACcCT(1)2#0.1 
FACT(2)2#0.1 
FACT(3!1=0 1 
FACT'(@)=0 1 
IF: 1COUNT.L 
DEL=c os 
FACTI‘1)=20 01 
FACT(2)=0 01 
FACT(3)=0 0° 
FACT(4)=0 01 

201 CONTINUE 
DELTA: 1)£4BS'CUESS(1,1)* 


Tetor GO TO 2 


4 


ESS FOR INITIAL COSTA 
a 


HOPT, LUZ 


HH, HMe ABSER,RELER, 1 


& 
R THE FIRST GUESS 


hip 


CTION ATTAINED’! 


ATIONS OF GUESS 


01 


DEL) 


ME VEC tO R oa 


4a) 
melee oe, 
AX Ye ea DIE SOIRIED = VEO) Si 
eae a) Pel 
S FOR COSTEATE VECTORS /a@4 2x E10 31) 


FLGLG! 


DELTA(2) 
CELTAI3) 
DELTA(4) 


=LBS(GUESS(1,2)#DEL! 
fABS'GUESS(1,3'*DEL) 
=ABS(GUESS(1,4)3DEL) 
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7 
°.eeeprety awed ae eee 


a a 
viens saneiaieless ae ne 
* .T tase wart 
-f 6h 26.@)) 400 7bd CeT RED CO toner ; fidi-y 


tee ere ae adi 
y 


- 


ero’) CO1tm ateen cme oe 46,90! #148, ¢ 


aa 2aaet. 


: “—™ 
comer pear? wees as hatad~ ea 
Bee ta 


aan 


no 


START ITERATIONS OF NEWTON-RAPHSON METHOD 


14 
WA 
Nes 


po 10 1=2,5 
Do We Jean, 4 
GUESS(],J)]8GUESS(1,J) 


GUESS(2,1"SGUESS(2,1)+DELTAI 
GUESS:3,2)=GUESS(32,2'*DELTA' 
GUESS'4,3°=GUESS(4,3'*DELTA\ 
GUESS(5,4;=GUESS(5S ,4;+DELTA! 


PUN 


BiG e278 
x20 

Vary eet TO 
Y(2)20. 

VP Veer 
Y(4&)=20 
yY(S);=Guess(1,i 
Y(6)8GUESS(1,2) 
Y(7)12GUESS(1,3) 
Y(&8)=GUESS(‘1,4) 
pO 12 K=1,MFIN 
BON V4 Lew 2 
DA=+5 

HH=+5 

HMXk=S 

RELER=0 1 
ABSER=0 

CHUL ROOTHY, HOPT, TOPT _DROPT Luz) 
TEMP=TOPT 

HAM=HOPT 

DHAM=DERHITEMP YY) 


CALL RKF(®X,Y,NOIM,FCT,DE, HR, HMX,ABSER ,RELER 


CONTINUE 

CONTINUE 

CALE OWTP TX, VY, PRPLAG, LU) 
WRITE: LUI, 15 °%,(Y OM! , Me, Bs 


15 FORMAT V//1X.F7.0,36,80/E10.3, 1x1 //) 
Path, ara C2 1) 
Par en=. 030) 
Fae, a ve Neca 
Se. 4 eyes") 
11 CONTINUE 
COMPUTE THE JACOBIAN 
BiO> VG as 4 
GO 6: =" 4 
16 QaGui wiv e cur, Pee ge Py ey DET Ard" 
WRITE CL UT esse dab ne oi, Wet So ee, 4 
33 FORMAT!4(/1x%,41£10.3,2Xx" 


INVERSE THE JSCOBIAN 


‘ec! 


ze 


24 


EiOl tia) Nee 4 

po 17 v=t,4 

MeTeCue 126 

AMATIMIZEJALTII,J 

NROW=6 

EPS=1 Ob-0O3 

CALL INVZ'AME&T,NROW, NROW, 4: FRI EPs Bik 
1L=46 

WRITE(LU1,28 1ER 


FORMAT(1X. “ERROR CODE OUTPUT OF INV2 1x 
WRITE COUT , sa vt AT ty 127 je 
FORMAT(61/1X%,4(010.3,2%°' 


C NEWTON-RAPHSON ITERATION 


(Star Men 


18 


30 


124 


FAC=1 D+00 


1S THE RELAXATION FACTOR 
Doe vs We, 4 
1H hfe ae rg al 
wooL=t 


cait MULTI(£1,NROW,NROW,BB NCOL NROW, FAC, 
WRITEULUN, v9" TER 

FORMET: 1X, 
Diores s)he, 
ctela= CG sek 
WRITE(LU1. 
FORMAT(//1 
Dio! 12:9) Jie 1 

GUESS, 17 


eFAC ht Ta 
Se Get 1) 
41010 


npeW- Db 
tw 
we 
“Nn 


Ew CUESS 
° 


m 
+z 


UESS(1 
UESS'‘1 
DiGi srs! 
UESS «1 


<<<«<«<<~< 
prone 
onnan000 


WRITE+LU1,3 
FORMAT(//1X, 1 TERATION NUMBER 
PNCENGy pay) o 


CuEIEOUNT 0 Mute, Les, 80) 


ICOUNT=1]COUNT?1 
Kz0 

Dat, 

Dior 23" 7 
po 124 L 
DA=+5 
HH=+5 
HMx=5 
RELER=O0.01 

ABSER=0 

CALL RDOTH(Y,HOPT,TOPT,OHOPT,LU2) 
TEMP=TOPT 

HAM=HOPT 

DHAM=DERH(TEMP,Y! 


,MFIN 
5 We 


nou 


“ERROR CODE OUTPUT OF MULTI’, 1X, 


~1FLAG: 


ey dik THE REW GUESS 


cat RKE IK. NDIM,FCT, DA, HH, HMX,ABSER RELERIFLAC) 


CONTINUE 

MM=K-L 

Lek 

CALL OUTPI(X,Y,I1FLAG,LU1! 

CONTINUE 

WRITEI(LU1,1251K,(V(d!',Jz1, 8) 

FORMAT‘1K, ‘RESULTS FOR THE NEW GUESS'/1X, 
EDOM NB EWON Sake 7) 


Coys eC at to gil rd cea 
U2 ag Wie 4a Ye 3b ir Jaa i 


IER! 


Sey 


/1x 


@ afe 
lene? 


anannn 


n 


Naanan 


6S 


6é& 


40 


500 


FH VFS es hs = 4 
Fiir,@izyvyis) 

BES AB SiR ta. Wel i: 

Ee eG Or02:S 7 ClO) TO) 6:8 
PRS abSu ety 20 

Tete Bee Ga O1O2'Sia SiO T:0) 6:8 
FFZABS(F(1,3)) 

HEGRE. GTS O. O255 GOs Ton 68 
FFSABS‘F:'1,4), 


TEMER GETS OC O2 Si COs TO 68 

WRITE(LUI,69) 

FORMAT‘ 1X, OBJECTIVE FUNCTION ATTAINED’) 
Geo TO 40 

CONTINUE 

co TO 1000 


CONTINUE 
WRITE(LU1,500) 


FORMAT(1X,’ITERATION HAS BEEN TERMINATED SUCCESSFULLY’) 


STOP 
END 


SUBROUTINE FCT(X,Y,DERY) 

DIMENSION Y:‘&:,DERY‘&;,YD(4 
RESL MO,K2,KP,K20,K1,KD,10 
COMMOh TEMP ,MO,10,KHOM,DHAM. YD 


COMPUTE THE DERIVATIVES 


2 


1 


3 
ao 


K20=100 &@8EXKP:-2960./TEMP;: 

K1=1.28E+*OS92EXP1-S962S /TEMP 
= 6 

Pe Tee 2 eR 2 Os 2 fF 

$6S0C /TEMP+250 

$000 /TEMP-185 

25, /7TEMP-C. 25 

= 

K 


" 
non 


COYOC2S P2YDI( 2) Ss 314+ BsV(2I2V(21FYD( 212 YDI(2)4€ 
ZOZ EXPER: 
2.8F2zkKD210 

Dee KO2 ¥ C2} 
r=SORT(K23S38Y 
PS2 Fee You is Oe YiD U3 4 
SQEKZEC(MOZTU-V¥O21FYD( 2p rez20/¥O1 4, 


Cea si A ay) Quits yi Pst 21 Ye pit > 
2 
1 
PSK o ys peSORTUK22 Sev TN 2 = iT =V 20s V Ores 1 /¥ D4 
2 
2 


=DERY/5'-S§2y 
PAZYU(2Z2V2ZVI2)27YD( 2)" YD(218YD(21492 sBzVYi2isyD(2) 
Sas + OK 2 haps ¥s Pz Die Ta LY. De 24 
=DERY‘6+13SORT:1S3)V¢6 118K 2:4Y 

PP ee (25) S¥ O02) + 2 sD Kis 

K2*?=MOSMO77' 8: /YD' 4) 

OERY:'6++DKb 


SN AQ Oil ge sree Gp ava 


ie ee 
monn om 


© 
° 


FTHgwM FE enn nnpuNn — 


SUBROUTINE OUTP‘£.Y¥,1}FLAC,LUI: 
DIMENSION 16) ,YDi6) 

RELL mMC,10 

COMMON TEMP ,.MC 10,HAM,DHAM. YD 
X3=4/60 

WRITE(LU1, 2/1 FLAG 

FORMAT!‘1XK, 16) 
WRITECLUN,1)X3.¢¥011,721,8), TEMP 
BORMAT OX, F6g 1, 2X) 80E 10.22% 1, FE 1) 
AMN=1002Y(2)2MO/V13)/YD(3'2YL'2) 
AMW=10027Y(4)/Y'2)/MO2YD(41/YD'2 
PD=QMW/AMN 

WRITE: LU1,3°2MK,AaMW, PC, HAM, DHAM 
EQGMATUGX 2 NE2O, 33%), FS. 2,213, E903 
RETURN 

ENC 


COMPUTATION OF THE HAMILTONIAN 


THIS CODE INTEGRATES A SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL 
EQUATIONS BY RUNGE-KUTTA-FEHLBERG METHOD WITH AUTOMATIC ESTIMATION 


OF 


FUNCTION H(T,Y) 

BaMEN SONS Ye) | yD a} 

REAL KipK2,K20, 16 ,MO,KD 

COMMOK TEMP ,.MO,10,HamM,DHAM, YD 

20S tOor a= EXP =29'60. /T } 

K1=1,28@E+OS2EXPi-9628 ./T) 

A=-96500 /T+250 

B=76000 /T-185 

G22 5e/1 -0 . 25 

Feo.e€ 

KD=(K1222)/K20/2 /F 

REAPCCY(212YD(2) :223/4BrV(2)3V1( 2127 YD( 212 YO(2)4C 
K2Z2=K 202 EXP (R) 

S=2. *2F2#KD210 

H=-Y(S'2KOrY( 1) 
HEH+SORT(S*K24Y1 1/38 (1-V(208YDI12)'/YO12i2 716) 
HMeHeS2zVYite2eVYDO(312Y!17) 


2 


ea 1ST 


HEH+2 =K2EMO2MO2(1-Y(2)8YD(2) 1261 -YC(2ISYDI(2))/YDI4)2V18) 


RETURN 
END 


SUBROUTINE RKF(A,Y,N,F,00,H,HMX,ABSER,RELER,IFLAG) 


LOCAL ERROR AND STEP SIZE ADJUSTMENT 


REAL Y¥(8),YTEMP(20),TEMP(20' ,R‘ 20! 

REAL K1120) ,K2120) ,K3‘20 ,KG120),K5(20!,K6(20) 

REAL U,RELER,ABSER,HMK,A,B,DL,HMAX,H,HKEEP,ARG,RATIO 
U=1.0E-12 

IF (RELER LT.0 O£LO.OR.ABSER.LT.O O£O) GO TO 18 
IFCRELER+ABSER.EOQ ©. .O£O) GO TO 18 

IF(HMX.LE ©.0€0! GO TO 18 


B=A+Da 

ASV=h 

IF(ABS'(DA).LE.13. OEO2UsAMAX1(ABS(A!),ABS(B:)'CO To 18 
HMAX=AMINTI)HMX,ABS( DA)! 

IF‘ABS(H) LE.13.0E0#UsABS(A)' H=HMAXK 

1ADUS=0 


H=SIGN(AMIN1(ABS(H),HMAX) ,DA) 
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o- M0 a _ 


Pe oe 


Sh a privet Pree 


t pens 1@ 
soho wee @ Pag = ise eae 


& @te OrOve@Gomoge @> are svete a : 
eh ; 


Gs = 45a @-—GF° 6 Pe ae 
Came 


i iecteeeats eee 
7 
6 say 


_ : PO ; 7 
‘eee a adudeanait cabepaye ass ' 
er ee oe ne oiee 


pees ees cuca 


361 
362 
363 
zee 
365 
36€ 
3e7 
BEE 
3és 
370 
371 
372 
ars 
374 
a7 Ss 
37€E 
377 
378 
378 
380 
3681 
382 
363 
3e4 
385 
386 
387 
36é 
3es 
z90c 
ao 1 
392 
3933 
394 
ass 
396 
ao7 
3ge 
3s¢ 
400 
801 
602 
4o3 
aoc 
40s 
aoe 
407 
BOE 
gos 
£10 
a1) 
412 
413 
G16 
ais 
aie 
a17 
416 
41s 
420 
421 
a22 
623 
426 
425 
G26 
427 
a2e 
429 
430 
43) 
432 
432 
434 
435 
436 
637 
4a3e 
4395 
44o 
rye 
442 
443 
446 
eas 
BSE 
4aa7 
448 
@as 
aso 
451 
452 
453 
454 
ass 
4aSE 
as? 
ase 
4ass9 
abo 
46) 
462 
463 
46a 
aes 
Gece 
467 
SEE 
aes 
470 
471 
472 
472 
ays 
475 
476 
a77 
476 
47s 
480 


(1. 25£0+10. EO#U)isABS(H)} CoO To 4 
HKEEP =H 
140NS=1 
H=B-£ 
4 Cte Fea Ka 
Ss CONTINUE 
Dd’ 6° 1=1.,N 
é YTEMP(CTIEY( 3 +46 25EOrHrK1( 1) 
ARC=A+0 25E08k 
CALL F: ARG, YTEMP K2: 
OP 7) AEM oN 
7 YTEMP(T)FY61)+HE(K1(1)*(3.E0/32.E60)4#K2(1)219.£0/32 EO)) 
ARG=A+H2(3 £EO/8 EO! 
CALL F: ARC, YTEMP,K3) 
OC <8) ea, N 
& YTEMPUT )EYC 1) He (K101)2(1832.60/21987 EO:-K2(1)381 
¥4K3(1)2(729€ £0/2187 E0)) BONA TDA SCRE OREN SESS 
ARC=A+H2/12.E0/13 EO): 
COLL FY ARG, YTEMP Ka) 
DO §$ J=1,N 
s VTEMP( I] )EY(1)¢H2 (K1(J 31835 EO/21E EO)-8 EO#K2(1) 
=+K3(1)*8(3680 £EO/513 EO!-K411)2(845.£0/4108 EO)) 
ARG=A+H 
CALL F‘ARG,YTEMP KS) 
00 10 1=1,N 
10 YTEMP(TVEV(T )+He i eK 31112 (8. £0/27.E0)+2.£0#K21(1)-K3'(1)*(35468,.E0/ 
2565 £O'+KG11)2(185S9 E0/4104 EO)-K5(1:2(11 EO/40.E0); 
ARG=A+O0 SEOFH 
CALL FI(ARG,YTEMP,K6) 
BO’ te! FST, N 
TEMP( 1] )EK1(1)8(25 £EO/216 .EO0'4K3(1)211408 .EO/2565 Eo) 
24KO (11212187 EO/6104 EO)-0O. 2E02KSi 1) 
aot YTEMPC] SV (0 J )¢ee TEMP) 
BO 12 Bet, N 
12 ROTVSRTOIGDYSIEO EOC-K3SITN2 1128 £O/4275 EO \-K8(1izt'2197 EO/75240.E0: 
SaKST HSC EOtKEt] Ie 12 EO/SS £0) 
RALTIO=0 OEC 
DOM aia) Ves .N 
TREABSIR' 11) /t(RECER®2BS(YTEMP:1):+ABSER: 
hes RATIO=AMAX1(RATIC,TR 
WEIR SteOne Gili £0 G0) 7:0) 4:5) 
OO) 74 Ms 05 'N 
14 YUTPEYTEMP LT) 
A=her 
HDA Se EO; wel GG; FO? ‘Ve 
RETIO=AMAXK1:‘(RATIC,€ SS3E6E-04: 
ve RETIO=AMINI(RATIC 40SE EO 
H=O BEOtH/SORTISORTIR&ATIO?) 
TFCABS‘H: LE 13. .EO#UZABSIDA GO TO 19 
Rkroon OR nee Om GO eto 3 
TAODUS=0 
eye) sy foe a} 
VE 1FLAG=1 
HZHKEEP 
RETURN 
18 IFL&G=3 
RETURA 
1S 1FLAG=4 
RETURK 
END 
c 
FUNCTION DERH'T, ¥? 
Re aiha Ss SUBROUTINE COMPUTE THE JOERTVETIVE SF THE 
C HAMILTONIAN H,SUBJECT TO V=1/T 
[2 
DOIMENSTORN Y'6& YO.s 
REAG Ki aR 2 MOO, K 20, KP. KD 
COMMON TEMP,MO,10,HAM,OHAM, YC 
Wieteeny T 
EDR=16€298 
E2R=2S60 
Reet. 28k Oo s/E xP = 9:6)2:9) sy) 
K202100 4#EXPi-2960 /T) 
FeO .6 
KD=(KV8225/K20/F/2 
A1l=-86S500 
B1=75000 
Ci=2-25 
A=h1/T*250 
B=B1/T-185 
Sst he 6: 2:5: 
REHEAT (Y(CZIFYD(2) B8BI+BeV(2stV( 2 FVD! 21/FYD(214C 
K2ZEK2Z2O0TEXPIR' 
S=22FeKkDz10 
DY1=-KDzY11) 
DYZZESORT(K29S2Y( V1 ;3161-V(212YD(2:+:/YD'2 
BYS= S25 (RV AYD C31 
DYO=2.8K2EMOZMOE!1-VYIZIFYOD(Z Bt 1-VYI2ZIFYDI2)/5/YO'4 
Fr:S-EOR*DY1 
F221 -E2R4A1F( (V1 2)8YD(2) )s2BIi+Brey( 2IFVI2IFYDIZ2IFYD(2) 
THOT -EDRI/2 sDYV2 
F3=-EOR#DY3 
FOL(-E2R+AIF((Y( 2128 YD(2) ) 883 eBI2ViZIFV(ZIVYD( 217 VDI 214011804 
DERHEFI2V(S)+F28V(61+F 3 V(71+F Oe V (8! 
DERH=-DERH2 VV 
RETURN 
END 
c 


TFC ABS(B-A? GT 


SUBROUTINE ROOTH(Y,HOPT, TOP 


C THIS SUBROUTINE FINDS THE ROOT 


C THE METHOD USED 


NnnNnNnanaANnn 


1S) AN HALF NTE 
DIMENSION Y(8) 
712313 
T2=363 
FISDERH(T1,Y) 
F2=DERH(T2,Y? 
P=FisF2 
El util th lib ee Sef 
32 CONTINUE 
3 WRITECE U2 4) 
4 FORMAT(1xK,’ERROR OCCURED IN 
1,’‘DH/OU HAS 2 ROOTS OR NONE 
[Mo saMt Vives 2a \ Ae Kod 
TESIS+FLOATII) 
HT=HIT,Y? 
DH=DERHIT,Y) 
9 WRITE(LUZ2,10)7,HT,DH,Y(2) 
10 FORMSLT(1X,F5.0,2X%,E10.3,2x, 
TEAM Tay S508 


jp leMaltell 50 LASi7 3a 
OF THE EQUATION DH/DU=0 
RVAL METHOD 


SUBROUTINE ROOTH’/1X 
IN THE INTERVAL (313,363'') 


EGO dis, FSi.) 


eS 


in seaman euler kh oe wnie 
digenes i er 


: ee i f 
439 ; 


~~ 


Ms i 
ne" Stew tt 
oa hee :, 
ore be # 
a 


Ss 


een ii. 
s1-00 pe ae pers 


- 7 


a 


2 


a" 
set 


ee anata mt suns as 


se d-qeni@ ges ssh eqeebmergens a. 
—_ th 
aria) 
oot) |) WeGe nt Fae 1 e lade 
Pe eee 0 doves wlan anenee 


 hctar aiialia iat of ae , oe 


a5) 7 TOPT=T2 149 


4E2 Gc To 40 

ae2 © TOP Tet 

ape GO TO 40 

aes 2 FPR ie fe. s 

EE 6 TOPTESET1 

4687 Go TO 35 

“ee 5 TOPT=T2 

aes Go TC 39 

480 1 CONTINUE 

ast c 

asi 200 PIs CTy+ T2772 

ag2 C TEST THE MAGNITUDE OF THE INTERVAL Ovanun ties 
age ERROR=ABS:7T2-7T1) 

49s IF (ERROR .GT.©.0001) GO TO 100 
aSE TOPTsT3 

ac7 GO To 39 

ass 106 F3I=DERH'T3,Y) 

ass F2=DERH(T2,Y) 

500 PHF3sF2 

501 HEPA ee es 

502 12 Tera S Ve, 1s 

503 16 TOPT=T3 

504 cc TO 39 

505 TS TOP b= TZ 

SOE GO Te 39s 

507 TY Ther? 

soe T2222 

508 GO TO 200 

510 Us Tee ty 

51) L=aT 3 

512 GO TO 20c 

§13 c 

s14 39 HOPTEH(TOPT, Y) 

515 OOFT=DERH'TOPT YY: 

S16 c WRIRECEU2, SEMTOPT HOP T DOP, yt2 
et 7 c So eORMrentOPK eee SUS kh , EO ai) 
sié &8O RETURN 

§1s END 


mor FILE 


¥ 
is 


+¢ 


7 
——> ane 
ped 


a 


; 


Fou, @re@_* as 


if 


a bears teks rr 


‘ 
- 


J 


3s 
> 


— 


wrmionuw run — 


10 


nanan 


no 


c 
c 
c 


nNnnn 


fe; 


DISCRETE CONTROL ALGORITHM 


EXTERNAL RKF.FCTY,FCTP,OUTPY,OUTPP,H, DERH 


DIMENSION Y(4),P\4),YD16) 


DCIMENSTON CONTR: 400) ,AUX1:14,400,12),6Hi 600) 


REAL MO,10 


COMMCN CONTR, AUKY,YD,MO,10,K,MFIN MSTP,¥ 


MSTP=12 
ALPHA=2 E+05 
10=0 01695 
Mo=9 3 
WRITE(6,980) 
8$O FORMAT/1X,°GIVE FINAL TIME’) 
READ‘ E,S3)MFIN 
83 FORMAT(13) 
TFEIN=FLOATIMFIN?: 
WRITE(6,91) 


$1 FORMAT(1X,‘GIVE DESIRED CONVERSJON,DESIRED MOMENTS’: 


READ(E,84i(YDI1),12=2,4) 
94 FORMAT(3F10.3) 


SELECT A PIECEWISE CONTROL 
WRITE16, 25: 


S85 FORMAT!I11X,’GIVE INITIAL ISOTHERMAL GUESS’) 


READ(6,96)TINI 
S86 FORMAT(FI1¢C 3: 
Oop 7 T=). MF TIN 
7 ECONTR CLT KET IN} 


WRT e 6: oy tateve Duet =z 
SRORMAR Uwe LK Vie DESTRED (= 
Seas DE STRED = 9 45.64 
I1COUNT=0 
1000 CONTINUE 
WRITE(6 ,&)1COUNT 
6 FORMAT(/1xX,’1TERATION NUMBER 


‘STE Pa 


PSU eG ae 
OF SGARY 2.679 


DESIRED = 


FeiCiya4 


INTEGRATE FORWARDS STATE EQUATIONS 


x=0 
yin a) 
Vv bi2 9 
Ysa 
yié) 
NDIM 
Lio 
DOs) kK 
Doe 12) a 
DA=+5 
HH=<+5 
HMx=5 
RELER=0 0} 
ABSER=0 
ChEL RKFtx,Y,_NOIM FCTY,DL,HH 
Pile 1) 3) aA 4 
Bie YC , 
Vz CONT VNWE 
MEK -L 
TEMUME DT ai) Gor To, 44 
LK 
CEL (OU TIPY UK GY 4 FEAG » 
11 CONTINUE 
WRT Ed6% '9°7 at ¥¢ 
WROD TET, 87 10 Yu 
S37 FORMAT‘ 1X ,4(E 


wun 
po0o0Oo- 


Sas 
Ge 

wun 
ie 
- = 5 


HM x 


ABSER RELEF 


STEP 3 
OBTALIN PITF) 
Pi1rizso 
Pi Sie sy 2 Na2) = Hy 0 
Pi Sis 2) SehiVd 39) 24 ) 
oli lke pg ee eS a ie ee | ‘ 
INTEGRATE BACKWARDS COSTATE EQUATIONS 
L=o 
Dakient 


K=TFIN260 

Dile 2 Hie Kea) ME oN 
DiOe (226 J= 1 MS THe 
DA=-5 

HH=-S5 

HMX=5 

RELER=0 01 
DBSER=0 


CALL RKF(X,P,NDOIM,FCTP,04,HH, HMX 


22 CONTINUE 
COMPLTE GRADIENT 
KKEMFIN-K+1 
MME=K-LL 
1F(MM.NE.0? GO TO 24 
GH(KK)=DERHIP) 
Lib = bib 
GO TO 25 
24 KKKEMFIN-LL*3 
GH( KK )=GHIKKK) 
25 CONTINUE 
M=K-L 
FLOM bet as BIOY TOs 2) 
CALL OUTPP(K,P,IFLAG: 
Lek 
21 CONTINUE 
WRITE'7,23)(GH'I),121 
23 FORMAT(/1K,5(E10.3,2X) 


Se Par6 
COMPUTE WORM OF GRADIENT 
GNORM=0 
DIGHHSis) DieciM iF dN 
31 GNORM=GNORM*(GH(])222) 
GNORM=SORT:‘GNORM! 
WRITE(€,33)GNORM 
33 FORMAT(/1xK, ‘’GNORM=°,E10 3/') 


FNORM=1AUXY(2,MFIN,MSTP)- 1 yes2+(AUXY(3,MFIN,MSTP!-1 


FNORM=FNORM+ (AUKXY(&,MFIN,MSTP)-1 


,ABSER,RELER, 


O);s22 


FSM Ue | 


IFLAG) 


I1FLAG) 


1232 


148 


121 
22 
ila. 
124 
125 
126 
127 
126 
12¢ 
130 
%31 
132 
133 
134 
135 
136 
137 
ase 
138 
140 
141 
142 
443 
144 
145 
146 
147 
14é 
149 
150 
151 
752 
453 
154 
155 
1SE 
157 
1seé 
155 
160 
161 
162 
163 
164 
165 
16E 
167 
16é 
tes 
170 
71 
172 
173 
174 
175 
V7E 
a77 
1768 
175 
180 
161 
12872 
163 
184 
yes 
186 
187 
188 
ree 
190 
181 
192 
183 
184 
195 
196 
197 
186 
iss 
200 
201 
2072 
203 
204 
205 
206 
207 
2068 
209 
210 
a] 
212 
mts 
214 
215 
216 
24° 
218 
21s 
220 
ee 1 
222 
223 
224 
225 
226 
227 
228 
228 
230 
231 
232 
233 
236 
235 
z36 
237 
238 
23° 
240 


WRITE: 6,34) FNORM 
36 FORMLT‘ 1x, ’FNORM=’, E10 3/1 
E1=ABS(Y(2)-1.) 
1FLE1.67T.0.05) GO To 101 
E2=ABS(Y(3'-1.) 
1FCE2 GT.©0.051.G0 te 101 
EXJ=ABSivigs-1 5 
HFhESWGT CC .O5% GoTo 104 
WRITE(5,103)(CONTR(I1),}=1,MFIN} 
103 FORMAT! 1X, 10F7.2) 
WRITE(E, 102) 


102 FORMAT(1X, OBJECTIVE FUNCTION ATTAINED’) 


CO TO 40 
101 CONTINUE 
C OTHERWISE COMPUTE & NEW GUESS 
DO 32 1=1,MFIN 
CONTR‘ T)=CONTRI I) -ALPHA2GH(1) 
32 CONTINUE 
WRITE(S,98)( CONTR‘ 1), 1 
9&6 FORMAT(1X,10F7.2) 


Aro NED 


c 
TCOUNT=1COUNT+1 
GO TO 1000 

c 

c 

40 CONTINUE 
WRITE(&,99)(CONTR(I?),}=1,MFIN« 
99 FORMAT LUX, 10F 7 . 2 

stop 
END 

c 

c 

c 

c 
SUBROUTINE OUTPY(x,Y,1FLAG) 
DIMENSICN Y'4),CONTR‘ 200: ,AUXY'14, 400,12) YD: 
REAL MO,1C 
COMMON CONTR, AUXY,YD.MO,1C,K MP UN UMS TP) oD 
X1=xK/60 
AMN=10C FY! 2)8YD(2)8MO/Y¥13//YD:i2 
AMW=100 fVi Ss /Y¥(2)/YD(2+/MO4rvDias 
PO=LMW/ DMN 
WRITETT VEXT, CYT), VE1,4),AMN,AMW, PD, CONTRIK) , JFLAG 

Pe ORME MAE Ce tere hy OWE tO a p2k ty Poe ey 2 KEG ew 5 OR eds 

RETURN 
END 

c 

c 
SUBROUTINE OUTPP‘X,P,IJFLAG! 
DIMENSION P44) ,CONTR' GOO: ,AUXY'4, 400,12: ,¥DiI45 
RELL 10,M0 
COMMON CONTR ,AUXY YD.MO,10,K,MFIN MSTP, J 
X1=¥/60 
KK EMFIN-K+1 
TEMP=CONTP IKK! 
HISHI TEMP, PY) 
WR E te 0a Xt Pa DN Pod, 4) RW oP LAS 

Dee ORM arm OE) ta ah 2 x EO (3 ip ae wl eae be WRAY Ve fel 

RETURA 
END 

c 

c 

c 
SUBRGUTINE FETY (X.Y OER Y 4 
DIMENSTON Y¥UAT ODERVY CS) YO(4) CONTRI A000) AUXY 4 , SOC , 12's 
REAL mMOC,K2,K20,K1,K0,1¢ 
COMMON CONTR,AUXY,YD,MO,10,K,MFIN,.MSTP,J 
TEMP=CONTRIK) 

c 


f EOMPUTE THE DERIVATIVES 
K20=100 4O2EXP1-2960 /TEMP: 
Kity ,2BESOSSEKPL-9629  /TEMP) 
=O 6 
KD=(K1222)/K20/2 /F 
=-$6500 /TEMP+25C 
B=7500C /TEMP-185 
Sra 2 Sy eee = Cm. 
K2=EK 202 EXP (A®((Y(2'18YO0(2)12323)4B" 
$=2 2F2kD210 


Y(2)"FV¥'2197YD(2)2YD(2)+C" 


DERY(112=-KD#Y011 
DERVESTSSORT (Ka22Se YC Prix tite YU2hs Y Dit 7 i/ Vota 
DERYGarsse yt TU/'Y Dt s4 
DERYCAVEQ2T EK 221 (MOF ISYIZVZYDI 21 12zZ27/YDt4s 
RETURN 
ENO 
c 
c 
SUBROUTINE FPETPCX) P, DER P » 
DIMENSION PI4:) DERP(4!,YD: 4), CONTR! 4800; AUXY(4,400,121,Y'8 
REAL MO, 16. K2;K20,K0,K1 
COMMON CONTR,AUXY,YD,MO,10,K,MFIN,MSTP, J 
KKEMFIN-K41 
JSJEMSTP-J4+1 
TEMP=CONTRIKK) 
DiOre yea ari 4: 
T1l=1+4 
YC T)SAUXY(I,KK,JJ) 
AVES, re Poy 
c 
c 
C COMPUTE THE DERIVATIVES 


K1=1.28E+*O097EXP:'-9629 /TEMP! 
K20=100 40#EXKP:1-2960 /TEMP) 
Feo 6 
KD=(K1222)/K20/2./F 
=-96500./TEMP+250 
B=7S000 /TEMP-185 
=-25 /TEMP-©.25 
K2EK208EXP(AR((Y(2)8Y0(2))223'+B 
S=2 *F2kKD210 
DERP(1)=KO7Y(S8) 
DERP(1)=DERP(1)-SORTIK22?S/Y(1)) ( 
DERP(1)=2DERP(1)°S#V¥1(7!/YD(3) 
DK=-(3.2A2VI218V(2)21(YD'(21323)42 


Vite eV he oe y Dil 2h sy ONG 2 Die iGo 


jeVi2 2 YO) ys VN G2 PY Dki2 4 


eB2Y(21)2YD(2)2YD'2!)/2 


DERP(2)21 *DK211-°Y¥12'2VD(21//YD(2) 


DERP12)=DERP(2'2SORT‘S2#Y(1)2K2'2Y(E) 
DKKESeYD(2+2(1 -¥(2)FYD(2))¢62DKE(1-¥(212YD(2))801 


DKK=DKK*#K2*MO8MO#Y(8!1/YD(4) 
DERP(2)=DERP(2)+0KK 
DERP(3)=0 


-¥(2)8YD(27)! 
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one . 
“er te ts 
7 ; “S 


“ai 
lass 


or oe wee arene. 
- : 0 ie aha 20 


nat tats aoe a 


1 

‘ 7 

c 4) set Peemes eer atapa,) tren: 
@ tem oi Wy Het " 


eee ee a 
‘pawn @ re ; 
\4-gttl een: “re A 


i: Gar : Pee 
He OR RE pee 


241 
242 
243 
246 
245 
246 
247 
24E 

24¢ 
250 
251 

252 
253 
254 
255 
256 
257 
256 
259 
260 
26) 

262 
263 
264 
265 
266 
267 
266 
269 
270 
271 

272 
202 
276 
ZzITS 
27€ 
277 
2768 
27S 
280 
281 
2&2 
2&3 
284 
285 
26E 
267 
28& 
28° 
280 
291 
222 
283 
2964 
2e5 
29€ 
287 
256 
299 
300 
301 
302 
303 
304 
305 
306 
307 
306 
309 
310 
at. 
312 
a3 
314 
315 
316 
a7 
316 
319 
320 
aa} 

maz 
323 
324 
325 
326 
227 
32 
32° 
330 
331 
332 
333 
334 
335 
336 
337 
3368 
335 
340 
341 
342 
343 
344 
345 
346 
347 
346 
349 
350 
351 

352 
353 
354 
oss 
356 

357 
358 
359 
360 


n 


nanannHn 


10 


12 


13 


14 


is 


16 


18 


1s 


DERP: 4'=0 
RETURN 
END 


COMPUTATION OF THE HAMILTONIAN 


FUNCTION H‘TEME PB) 


OIMENSION Pi4s,¥(8@),¥YD(&), CONTR‘ 400: AUXKY\ 4, 400 12) 


REBL K1,K2,K20,KD,10,mM0 


COMMON CONTR, AUKY,YD.MC,10,"% .MFIN,MSTE, J 


KKEMFIN-K +1 


JJtt 
DOeis hz? ,4 
T1=1+4 


Y(IT)SAUKY(I,KK, Jd) 
Bk (1 aT Ba 
K20=100 402EXP\(-2960 /TEMP) 
K1=1.28E+*OS92EXP:1-9629 /TEMF) 
&=-96500 /TEMP+250 
B=75000. /TEMP-185 
25 eA eMe.- OL 25 
Fro 6 
KD=(K1®*2)/K20/2 /F 


K2=K2Z02EXP(ASY(2)9V(2)8Y(2)2(YD(2) a5 3 +Bevi2izyiz2is2yYDi2)2YD 2)4+D: 


S$=2.2FeKDz2]0 

O1=-KOzYi1) 
D2=SORT(S*K2=Y(1)/)8(4-V¥(2:2YDI2:1/YD(21 
D3=St#yYiir/YD'I3 


DO=2.tKZ2=MOFMOFI(1-VIZ2ISYDIZ2+/18(1-VI2+FYD'2/)/YD16) 


H=012Y(5)4022Y(6)4¢032Y(7)/4+DG2v/(8) 
RETURN 
END 


SUBROUTINE RKF(L,Y,N,F,DL,H,HMX ABSER RELEF,IFLAG) 


THIS CODE INTEGRATES A SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL 


EQUATIONS 6 


OF 


LOCAL ERROR AND STEP SIZE ADUUSTMENT 


READ yo wre Me i 2O0n TEMP C20 Rt? Ord 


RUNGE-KUTTA-FEHLBERC METHOD WITH AUTOMATIC ESTIMATION 


REAL K1(20; ,K2120),K3'20: ,K4(20:,K5!20),K6(20) 
REAL U,PELER,ABSER,HMX,4,B,0L,HMAX,H, HKEEP ARO. RATIO 


U=1 OE-O€ 


IFCRELEF LT CGC O£O OR. ABSER.LT © OEO) GO TO 18 


PETRECERtARBSER EQ © OFO)! GO TO 18 
TES HMA LEO LOGO} GO TO 18 


B=h+0L 

aSvee 

IFCABS( DA? LE.13 OFO2 UZ AMAX1(ABS( A+, ABS‘B'11GO TO 18 
HMAX=AMINT'HMX ABS: DL: : 

TFA‘ ABS(H'. LE 13 OEO*UZABS(A:' HE=HMEK 

TAODIS=0 


H=SIGN'‘AMIN11ABS(H HMAxX) DL) 
TFCABS1B-A}.GT (1.2560*10 EO™UIZABSIH)») 
HKEEP=H 


1ADUS=1 

H=B-L 

CANE Fata Ve Reatey: 

CONTINUE 

DO 6 11.N 

YHUEMP Yee Vi ly v+0" 2S E02 Msi tT 2 


ARG=L+0 Z25E03n 
CALL Fi ARC, YTEMP ,K2: 
OO 7 J=1,N 


GO TO 4 


Vine Mirena Olas Mane tits tat SFOs? OMe K24 hist 9) BOY S:2), EON 


ARG=A+H213 EO/E. EO! 
CACE FrARG, YTEMP , KS) 
DiGi es kai iN. 


Vem Tey Tis His OK tet pay es2, EO 2197 FOr 


SKS a tne 9 6 ie Or 2 1 8ir EO } 
ARG=A+H2(12.E0/13 EO! 
C&LL Fi‘ GRG, YTEMP,KS) 


Keg Meee Cn FOV rg 1/97) 21E 


BO| 9 t= 4), .N 

YTEMPC IT DSV'1 4H (K101121439 £EO/216 EFEO)-& EO2K2:} 

*K3I 11) 03680.60/513.FEO)-K40) 121885 .E0/41048 EO!) 

ARC=A+¢H 

CALL FrarG -¥Y TEMP. KS 5 

Bios VO t=), .N 

YTEMPUTDEYCT '4Hst-K1(1)218 £€0/27.£0)42 EO2K2'1)-K3i1)2:3584 EO/ 
2565 .£0:4K4(1)21185S.E0/4104 EO'-KSi(1)2111.£0/40 EO': 

ARG=A+0 SEO#H 

CALLE FrARG, YTEMP KE) 

OO it kent, 8 

TEMP11)=K1(1'2(25 £EO/216 £O'4K3(1'3(140& EO/2SE5 EO) 
oK4'1'212197 £O/8108.E0)-0.2E02KS (1) 

YTEMPi J] '=¥11)¢H= TEMPI 1) 

DO 12 J=1,.N 

R(1)2K111+1/360, £EO°K3(1121128 EO/4275 EO'-K4'112(2187.E60/75240 EO: 


*KS(1)/50 EO*KE(1)8#(2,E0/55.E0) 
RATIO=0.0EO 
COCs aN 


TREABS(R(1))/(RELER®4BS(YTEMP(1))+ABSER) 


RATIO=AMLX1/(R4ATIO, TR 
DFU AT POG Tat. FOU GO! 1:0 15 
DO 146 12=1,WN 

YCl eYTEMPI I) 

B=A+H 

IF: TADIJS.EO.1) GO TO 16 
RATIO=AMAXK1(RATIO,6 SS3IEE-04) 
RATIO=AMIN1(RATIO,4096 EO! 
H=O, BEO2H/SORT(SORT(RAT1IO)! 
IF(ABS(H) LE.13.£0*U7ABS‘DA!) GO TO 19 
IF(RATIO.LE.1.£0) GO TO 3 
1ADJUS=0 

Go Tc § 

IFLAG=1 

H=HKEEP 

RETURN 

IFLAG=3 

RETURN 

IFLAG=4 

RETURN 

ENC 


FUNCTION DERH(P) 
REGL K1,K2,K20,KD0,M0,10 


CIMENSION CONTR(400;,AUXY(4,400,12),P(41,Y014),Y(8) 
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os at 62) 0's aieas tr anekin 
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-@ dprie Wee cnay wiooneeet : 
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- - >. beru re 


-. rn onl ti Raye oP meee. 


6 ee 
“ “vieies oooh me pepe ab Fela at sibise o 4. j 
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CPRtagr oR Me eR 
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cs bg ret 
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361 COMMON CONTR,AUKY YD,MO,10,h ,MFIN 
362 KKEMFIN-K+14 ' PMSITP hc) 151 


362 anrsiee5) 

366 TEMP=CONTRIKK?) 

365 po 1 J=1,4 

366 ViT'SAUKY(Y ,KK,JJd! 

367 Tl=1¢4 

368 ve sok ee PHT 

36S C THIS SUBROUTINE COMPLTE THE DERIVATIVE OF THE 

370 C HSMILTONIAN &K,SUBJECT TO TEMPERLTURE 

37) c 

372 EDR=16298 

373 E2R=2960 

376 K1=1.28E+O092EXP(-9629 /TEMP) 

375 K20=100.407EXP(-2860 /TEMP) 

276 Feo.6 

377 KD=(K1882)/K20/2./F 

376 A1=-96500. 

3798 Bi=+7S5000 

380 Di=-25 

361 h=2h1/TEMP+250 

382 B=B1/TEMP-185 

3283 D=D1/TEMP-0.25 

384 K2=K2027 EXP. as = 2 

Ee oe Doe nr ce YC2IFV(2)8Y(2)2(YD(2)88314B2Y12)29(2)8YD(2)2YD(214D) 
386 OvY1=-KDzY11) 

we? DYZ=SORT(K27S4Y( 11) 8(1-¥(213YD(2))/YD(2) 

35é& DY3=S*vY(1)/YD(3) 

36° DY4=2.*K2*=MO*MOFI1-Y(2)2 Cua . > 
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91 FQ Wy 22 k= + 21 DP 
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333 F3=-EDR2DY3 

394 F4='-E2R4 A126 (VY(2)2YD(211223 1/4818 

395 Re ae et ese cca iter Ot ale Rpt e ye 
BSE DERR=-DERH/ TEMP /TEMP 

397 RETURN 
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EXperimental equipment was built in the polymer labor- 
aeory.or the sDepartment of Chemical Engineering at the Univ= 
ersity of Alberta, in order to verify the optimal temperature 
policies for a polymerization reaction, and particularly for 
the polymerization of methyl-methacrylate, as part of another 
project. The experimental apparatus is described in Figure 
Eines 
Te) Measurements. In order to verify the validity of the 
mathematical model and the optimal temperature policies, 
converison and molecular weight distribution must be measured 
Continuously during the reaction. Samples are therefore 
Pakenein a continuous manner from the bottom of the reacron 
by means of displacement pumps then directed through a 
density meter and a gel permeation chromatograph (GPC), 
Which are interraced to a HP 1000 computer. feasurements 
of conversion and molecular weight distribution are recorded 
Gpecharts ana in a data file of the HP 1000. 
Ti, Temperature Control of the Reactor. Since polymerization 
is highly sensitive to variations of temperature in the 
reactor, it is important that there be good CONtLOlen ene 
temperature in the polymerization reactor. In this experi- 
mental system, the temperature Vee control leds by “age Ly 
temperature controller. The manipulated variable of the 
temperature loop is the cold water flowrate entering the 
cooling coil in the reactor. Theuset) DOlnts OL snes TD 
controller can be activated manually or by means of a 
supervisory control program on the HP 1000 computer to 


follow the optimal temperature policy. The high temperature 
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necessary to start the polymerization is obtained by means 
of a steady flow of hot water in the outside jacket of 
the reactor ee Regulation of the reactoretemperature ac.gooc 
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